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Capítulo 1 



Introducción 



El principal objetivo de este trabajo es presentar una adaptación de los métodos de 
volúmenes finitos utilizados en la resolución de problemas provenientes de los procesos de 
sedimentación de suspensiones floculadas (o sedimentación con compresión). Esta adap- 
tación está basada en la utilización de técnicas de multiresolución, originalmente ideadas 
para rebajar el costo computacional en la resolución numérica de leyes de conservación 
hiperbólicas, en conjunto con esquemas de alta resolución. 

Se introducirán los métodos utilizados para la resolución numérica de leyes de conser- 
vación y ecuaciones parabólicas y la importancia del algoritmo de multiresolución en la 
aplicación de estos métodos. 

Leyes de conservación hiperbólicas 

Los sistemas de leyes de conservación son modelos matemáticos para situaciones físicas 
en que la cantidad total de la variable no varía con respecto al tiempo. En este tipo 
de situaciones, la cantidad de una variable física contenida en una región acotada del 
espacio sólo puede variar debido al flujo de la variable a través de la frontera de dicha 
región. Esto puede traducirse en una formulación integral que, bajo ciertas hipótesis de 
regularidad, se convierte en un sistema de ecuaciones en derivadas parciales. Si se toma el 
caso unidimensional (en espacio), las ecuaciones correspondientes son de la forma 

d t u(x,t) + d x f(u(x,t)) = 0, (1.1) 

donde u : M x M — > M m es el vector de variables conservadas o variables de estado, y 
/ : M m — ► M m es el vector de flujos. En problemas de dinámica de fluidos, estas variables 
son densidad, momento y energía. 

La ecuación (1.1) está provista de condiciones iniciales y posiblemente condiciones de 
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frontera en el dominio espacial acotado. 

Un ejemplo clásico para ilustrar el comportamiento de las soluciones en leyes de con- 
servación, es el problema de Riemann en un tubo de shock: dinámica de los gases. Se tiene 
un tubo lleno con gas, inicialmente dividido en dos secciones por una membrana. El gas 
tiene densidad y presión, en reposo, más alta en una mitad del tubo que en la otra. En el 
tiempo t = se rompe la membrana y el gas ñuye. Si se supone que el flujo es uniforme 
a lo largo del tubo, la variación se produce sólo en una dirección y pueden aplicarse las 
ecuaciones de Euler en una dimensión. 

La estructura de la solución del problema de Riemann implica tres ondas distintas que 
separan regiones en las que las variables son constantes. La onda de choque se propaga 
hacia la región de más baja presión; a través de esta onda, la densidad y la presión asumen 
valores más altos y todas las variables son discontinuas. Luego aparece una discontinuidad 
de contacto, a través de la cual la densidad es discontinua, pero las demás variables son 
constantes. La tercera es la onda de rarefacción (recibe este nombre debido a que la 
densidad del gas decrece cuando esta onda pasa a través de él) que se mueve en dirección 
contraria a las otras dos y tiene una estructura diferente: todas las variables son continuas 
y presentan una suave transición [21]. 

Ecuaciones parabólicas 

Se quiere estudiar un problema de valores iniciales para una ecuación parabólica. Para 
ello, para (x,t) £ $7 x [0, oo[, considérese la ecuación 

d t u(x,t) + d x F(u(x,t),d x u(x,i)) = S(u), 
u(x, 0) = uo(x) 

donde ahora el flujo F incluye a la derivada de u y este se define por un operador diferencial 
con difusividad constante v > 0, es decir, 

F(u(x,t),d x u(x,t)) := f(u) - vd x u(x,t). 

Se tienen versiones lineales y no lineales. Para la ecuación de convección-difusión unidi- 
mensional, se tiene 

f(u) = cu, 
S(u) = 0, 

con c > 0. Este tipo de ecuaciones es de gran utilidad, por ejemplo, para calcular el 
transporte de sedimentos así como el transporte de constituyentes en estudios de calidad 
de agua [12]. 
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En el caso de la ecuación viscosa de Burgers unidimensional, se tiene 

/(«) 



v 2 



2 ' 

S(u) = 0, 

Esta ecuación es un modelo sencillo para la propagación de fluidos, tomando en cuenta 
que existe viscosidad constante en el fluido. 

Para la ecuación de reacción-difusión (a > 0,/3 > 0), 
/(«) = o, 

S(u) = — (1 — u) exp 



2 v ' ^ a{l-u) - 1 

Esta ecuación conduce al modelo unidimensional de la propagación de llama premezclada 
[32], donde las difusividades de masa y calor son iguales. La función u representa la 
temperatura adimensional, que varía entre y 1, y la masa parcial de gas sin quemar es 
representada por 1 — u. 



Ecuaciones parabólicas fuertemente degeneradas 

Considérese una ecuación parabólica de la forma 

d t u + d x f(u) = d 2 xx A( U ), (1.2) 

con (x,t) e]0,l[x[0,T[ y 

pu 

A{u) := / a(s)ds, a(u) ^ 0. 
Jo 

En general se permite que a(u) sea cero en incluso un intervalo [0, u c ], en el cual la ecuación 
es de naturaleza hiperbólica, y a(u) es discontinua en u = u c . Dada la forma degenerada 
de a(u) y la naturaleza generalmente no lineal de f(u), las soluciones de la ecuación son 
generalmente discontinuas y es necesario considerar soluciones entrópicas. 

Una ecuación de convección-difusión fuertemente degenerada, con una función de flujo 
no necesariamente convexa que depende del tiempo, asociada a ciertas condiciones iniciales 
y de frontera como (1.2) se considera como el modelo clásico para los procesos de sedimen- 
tación-consolidación. La sedimentación es, a grandes rasgos, un proceso en que partículas 
o agregados son separados de un líquido bajo la acción de la fuerza de gravedad. Este es 
probablemente el método industrial a gran escala más importante utilizado en química y 
minería [31]. En soluciones relativamente diluidas, las partículas no se comportan en forma 
discreta sino que tienden a agregarse unas a otras durante el proceso de sedimentación. 
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Conforme se produce la floculación, la masa de partículas va aumentando, y se deposita a 
mayor velocidad. La medida en que se desarrolle el fenómeno de ñoculación depende de la 
posibilidad de contacto entre las diferentes partículas, que a su vez es función de la carga 
de superficie, de la profundidad del tanque, del gradiente de velocidad del sistema, de la 
concentración de partículas y de los tamaños de las mismas. El efecto de estas variables 
sobre el proceso sólo puede determinarse mediante ensayos de sedimentación. Esto hace 
que sea de gran utilidad en la modelación de estos fenómenos, la teoría de problemas 
inversos (ver [2, 16] entre otros). 

Desde hace ya varios años se ha estado trabajado con mucho énfasis en mejorar los 
fundamentos de los modelos existentes para este tipo de procesos. Grandes avances se 
deben al trabajo de Bürger et al. [4, 5, 7, 8, 9, 10] entre otros. Para una descripción 
detallada de estos procesos y su modelación, se recomienda consultar [4, 11]. 

Por las características de este tipo de ecuaciones, no es posible aplicar ni la teoría de 
ecuaciones estrictamente parabólicas, ni la teoría establecida de soluciones de entropía de 
leyes de conservación [10]. 

Método de multiresolución: Motivación 

Generalmente, el vector de flujos en una ecuación hiperbólica o parabólica, está formado 
por funciones cuya dependencia de las variables de estado es no lineal y esto hace que no 
sea posible deducir soluciones exactas para estas ecuaciones. De aquí nace la necesidad de 
diseñar métodos numéricos que aproximen convenientemente estas soluciones. Este es un 
problema general que afecta a la mayor parte de las ecuaciones en derivadas parciales no 
lineales, sin embargo, existen razones para estudiar esta clase particular de sistemas: 

■ Muchos problemas prácticos en ingeniería y ciencia involucran cantidades que se 
conservan y conducen a problemas del tipo ley de conservación. 

■ Existen dificultades especiales y específicas a esta clase de sistemas (por ejemplo la 
formación de ondas de choque) que no se observan en otros problemas no lineales y 
que deben tenerse en cuenta en el diseño de métodos numéricos que aproximen sus 
soluciones. 

■ Aunque se conocen pocas soluciones exactas, la estructura matemática de las ecua- 
ciones y sus soluciones es cada día más estudiada. Este conocimiento se puede apro- 
vechar para desarrollar métodos adecuados a las características de estos sistemas y 
sus soluciones. 

El hecho de que las soluciones de este tipo de ecuaciones admitan discontinuidades plantea 
varios problemas, tanto desde el punto de vista matemático como numérico. Es evidente 
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que una solución discontinua no puede satisfacer la ecuación en derivadas parciales en el 
sentido clásico. La teoría de distribuciones provee de una herramienta matemática muy 
útil, pues permite caracterizar las discontinuidades admisibles y definir el concepto de 
solución débil de un problema diferencial. 

Sin embargo, la clase de funciones continuas a trozos es demasiado amplia para garan- 
tizar unicidad de solución. Generalmente existen soluciones débiles con los mismos datos 
iniciales. Puesto que estas ecuaciones son modelos para situaciones físicas reales (o al me- 
nos esa es la motivación), es obvio que sólo una de estas soluciones puede ser aceptable 
desde el punto de vista físico. El hecho de que existan otras soluciones espúreas es con- 
secuencia de que nuestras ecuaciones son tan sólo un modelo que ignora algunos efectos 
físicos, particularmente en el caso de leyes de conservación, los efectos difusivos y viscosos. 
Aunque estos efectos (y otros) pueden ignorarse en la mayor parte del fluido, cerca de las 
discontinuidades juegan un rol esencial. 

Estas consideraciones conducen a la imposición de determinados criterios basados en 
consideraciones físicas que permiten aislar la solución físicamente relevante entre todas 
las posibles soluciones débiles. Este tipo de criterios se conocen como condiciones de en- 
tropía de nuevo por analogía con la dinámica de gases (en este caso, la segunda ley de la 
Termodinámica: La entropía nunca decrece). En particular cuando las moléculas del gas 
pasan a través de una onda de choque, su entropía deberá aumentar, y esto proporciona 
el principio físico adecuado para determinar de manera unívoca la solución con sentido 
físico. 

La aproximación numérica de este tipo de soluciones incorpora un nuevo conjunto de 
problemas. Las discretizaciones de la ecuación en derivadas parciales mediante diferencias 
finitas ocasionarán problemas si las soluciones que se quieren aproximar son discontinuas. 
Estos problemas son de dos tipos. En general, los métodos numéricos de primer orden 
incorporan difusión numérica; esto facilita la convergencia a la solución entrópica, pero 
limita la utilidad real de estos métodos. Los métodos clásicos de orden superior reducen 
la viscosidad numérica pero incorporan términos dispersivos y dan lugar a oscilaciones 
numéricas que pueden desencadenar inestabilidades no lineales o hacer que las aproxima- 
ciones numéricas no converjan a la solución físicamente relevante. 

Los esquemas numéricos diseñados para la aproximación de las soluciones de este tipo 
de ecuaciones deben poder escribirse en forma conservativa. Esto garantiza que si las apro- 
ximaciones numéricas convergen, lo hacen a una solución débil de la ecuación (Teorema 
de Lax-Wendroff). 

Si un método conservativo satisface además algún análogo discreto de las condiciones de 
entropía, el límite de las aproximaciones numéricas será precisamente la solución relevante 
desde el punto de vista físico. 
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Una excelente clase de métodos conservativos para la aproximación numérica de las 
ecuaciones hiperbólicas y parabólicas, son los métodos de alto orden de precisión. Estos 
proporcionan perfiles bien delimitados y sin oscilaciones cerca de las discontinuidades. Un 
aspecto importante a tener en cuenta de los métodos de alto orden de precisión, es su 
elevado costo computacional, el cual es aún mayor bajo las siguientes condiciones: 

■ Sistemas de ecuaciones. 

■ Más de una dimensión. 

■ Un gran número de puntos en la malla. 

■ Extensos períodos de simulación. 

Método de multiresolución: Descripción 

El método de multiresolución es una técnica destinada (al menos, originalmente) a 
rebajar el costo computacional asociado a los métodos de alta resolución. En situaciones 
estándar, el comportamiento de la solución w(x, t) como función de x es altamente no 
uniforme, con fuertes variaciones en regiones puntuales y un comportamiento suave en la 
mayor parte del intervalo computacional. La técnica de multiresolución (al menos, en la 
forma en que será utilizada en este trabajo) fue diseñada originalmente por Harten [26] 
para ecuaciones hiperbólicas y utilizada por Bihari [3] y Roussel et al. [32] para ecuaciones 
parabólicas. Se desea estudiar la aplicación del método de multiresolución a los métodos 
existentes para modelar fenómenos de sedimentación de suspensiones floculadas [11]. 

Dado un método en forma conservativa y una malla uniforme apropiada para la solu- 
ción numérica del problema de valores iniciales para una ley de conservación hiperbólica 
escalar o una ecuación parabólica, el método de multiresolución aproxima la solución a una 
tolerancia prescrita de una forma más eficiente, entendiendo por eficiencia una reducción 
en el número de veces que se calcula el flujo numérico con el método de alta resolución. 
Para ello se consideran los valores puntuales o medias en celda de la solución numérica 
mediante un proceso jerárquico de mallas anidadas diádicas, en el cual la malla dada es 
la más fina, y se introduce una representación que contiene la misma información. 

La representación de multiresolución de la solución numérica está formada por sus 
valores puntuales en la malla más gruesa y el conjunto de errores por interpolar los valores 
puntuales de cada nivel de resolución a partir de los del nivel próximo más grueso. La 
compresión de datos es realizada haciendo cero las componentes de la representación que 
están por debajo de una tolerancia prescrita, e incluso eliminando de la malla a los puntos 
cuyos errores son menores a esta tolerancia prescrita; por consiguiente en lugar de calcular 
la evolución en tiempo de la solución numérica en la malla dada, se calcula la evolución 
de su representación de multiresolución comprimida. Como la transformación entre una 
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función y su representación de ondelette es rápida, la proposición de efectuar la gran parte 
de los cálculos en la representación de multiresolución es factible y atractiva. 

La información contenida en el análisis de multiresolución de la solución numérica es 
utilizada para identificar la localización de las discontinuidades en la solución numérica, 
y diseñar métodos que mejoren el cálculo del flujo numérico. Esta información es de gran 
utilidad al momento de calcular los flujos, pues el procedimiento correspondiente toma 
en cuenta la regularidad de la función. Además, la eficiencia computacional del método 
de multiresolución está directamente relacionada con la razón de compresión de los datos 
iniciales, es decir, la solución numérica en la malla más fina [ ] . La eficiencia del algoritmo 
se mide mediante la tasa de compresión y el tiempo de CPU. 

Programa 

Este trabajo se organiza del siguiente modo: En el capítulo 2 se revisarán los conceptos 
básicos necesarios para el análisis de multiresolución propuesto por Harten [26]. En el 
capítulo 3, se utiliza este análisis para desarrollar un método de alta resolución en mallas 
generadas por multiresolución, diseñado por Kozakevicius y Santos [29], el que será apli- 
cado a leyes de conservación hiperbólicas escalares. Se muestran resultados de los test 
numéricos realizados. En el capítulo 4 se analizan las ecuaciones parabólicas escalares y 
un método numérico que utiliza la multiresolución y la alta resolución (esquemas ENO de 
segundo orden y esquemas Runge-Kutta de segundo orden) como herramientas principales. 
Se utiliza una nueva estructura de datos desarrollada por Cohén et al. [1 í]. Se muestran 
los resultados de los experimentos numéricos realizados, coincidentes con los resultados 
obtenidos por Roussel et al. [32]. En el capítulo 5 se presentan los supuestos básicos para 
el problema de la sedimentación, analizando varios casos test. Se simula un proceso de 
sedimantación tipo Batch y se muestran resultados obtenidos aplicando métodos de mul- 
tiresolución a los esquemas desarrollados por Bürger et al. [5, 7, 8, 9, 10]. Se observa que 
el método de multiresolución es de gran ayuda para reducir el costo computacional en este 
tipo de problemas sin afectar la calidad de la solución. 
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Capítulo 2 



Multiresolución y compresión de 
datos 

En este capítulo se presentan los conceptos y definiciones básicas introducidas por 
Harten [26] para el análisis de multiresolución. Se presentan además herramientas adicio- 
nales utilizadas por Kozakevicius y Santos [29] para el desarrollo de métodos con mallas 
generadas mediante análisis de multiresolución. 

2.1. Análisis de multiresolución para valores puntuales 

Considerar Nq = 2 n ° valores 

«° = {«?}&' (2-1) 

correspondientes a los valores puntuales de una función u(x) sobre una partición uniforme 
de [-1,1]: 

G° = ,x Q ó = -l + j-h L , ^o = J-, «? = «(*?), Kj<N . (2.2) 

Se supone que u(x) es 2-periódica. Sus valores fuera de ]-l,l] son los de su extensión 
periódica: u® = u¡jy- , etc. 

Considerar el conjunto de mallas anidadas diádicas G k , k = 0, . . . , L: 

G* = {x}}£* , ,x$ = -l + j-h k , h k = 2 N *+ 1 h , N k = ^, (2.3) 

donde el nivel k = corresponde a la malla original, que es la más fina; y k = L corresponde 
a la malla más gruesa. Notar que G k está formada a partir de la malla más fina G fe_1 
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eliminando las componentes de la malla con índice impar, es decir 

'2j-1/j=1' x j - x 2j 



G k - l \G k = {x k ~\}% x) = x k 2j \ 0^j^N k . (2.4) 
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Figura 2.1: Diferentes escalas de valores puntuales 
Además se definen 

uj = u(a$) = u(x° 2kj ) = u° 2kj , 0^j^N k , (2.5) 
por lo tanto este proceso (ver figura 2.1) permite obtener u k a partir de u k ~ l mediante 

u) = u\~\ l^j^N k , (2.6) 
u k ~ l -u k = {u k j\}%. (2.7) 

Sea I(x, u k ) una función de interpolación de la malla fc-ésima, es decir, 

l(x k ,u k ) = u k , 0^j^N k , (2.8) 

que puede utilizarse para obtener aproximaciones para los valores ausentes en la malla 
k — lésima 

ñ k j\ = I{x k jl x , u k ) , < j < N k . (2.9) 

Sea D k (u°) = {Dj}^^ la sucesión de errores de interpolación al predecir los valores 
puntuales de cada nivel de resolución a partir del próximo nivel más grueso 

D k = v*j\ - ü k r\ = u k j\ - l{x k -^u k ), 1 < j < N k . (2.10) 

Estos D k se conocen como coeficientes de ondelette o detalles. Es sencillo comprobar que 
los conjuntos de datos (u k ,D k ) y u fc_1 contienen exactamente la misma información, 

u k ~ l «-► (u k ,D k ) (2.11) 
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en el sentido de que existe una transformación uno a uno entre ambos conjuntos (notar 
que la cardinalidad es la misma: N^-i = 2Nk). 



Claramente utilizando (2.11) sucesivamente para 1 ^ k ^ L, se obtiene 

vP^iu 1 ^ 1 ) <-» (D\(D 2 ,u 2 )) = (D\D 2 ,u 2 ) <-> ■■■ 
<-» (D\D 2 ,...,D L ,u L )=:(u M ) T 



(2.12) 



donde um = D 2 , . . . , D L , u L ) T es la representación de multiresolución de u°, equiva- 
lente a la representación original. Esta permite extraer información sobre la suavidad de 
la solución a partir de los errores de interpolación. La transformación uno a uno entre u° 
y um 



u M = Mu 1 



■u 



° = M- 1 u M 



(2.13) 



es lineal si I(-,u k ) es independiente de los datos. En principio, puede utilizarse cualquier 
técnica de interpolación para X. En este caso se utilizará interpolación central polinomial 



T(x,u k ) = qj(x), x <E Ij = [xj-!,Xj], j = 1, 



(2.14) 



donde qj(x) es un plinomio de grado r = 2s unívocamente determinado por los datos 
{iij_ s , ■ ■ ■ ,Uj +s _ 1 ) en los puntos (xj_ s , . . . ,x k +s _ 1 ); el valor en x^J^ se calcula a partir 
del polinomio de grado r — 1 (es decir, cada esténcil está formado por r puntos consecutivos 
de la malla) que interpola los puntos (uj_ s , ■ ■ ■ , Uj +S _ l ), por consiguiente 



?! 



2j-l 



X(x 2 j_i, u ) 



con 



2 
4 



Ver detalles en el apéndice A.l. 



X>(4m-i + 4-')' 



i=i 



01 = 1/2 

01 = 9/16, 02 



1/16 



2s, 



(2.15) 



(2.16) 



En este caso M es un operador lineal que puede ser representado por una matriz de 
iVo x A^o- Sin escribir la forma explícita de esta matriz, se sigue que um = Mu puede ser 
calculado mediante el siguiente Algoritmo de Codificación 



u M = Mu u 



FOR k 


= 1,2,.. 


.,L 




u) 


- u* -1 

- u 2j > 


1 < 3 < N k , 




D) 


- U 2j-1 


-Et=iA(«í+,-i + «}_,) 


i < i < ^fe 



(2.17) 
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y además u — M 1 um puede ser calculado mediante el siguiente Algoritmo de Decodifi- 
cación 

f FOR k = L,L 



u° = M~ l UM < 



uiy = u), 1 < j < N k , (2.18) 
«27-1 = £?=i A(«J+i-i + «J-j) + i < i < N k . 



Notar que el algoritmo de Codificación va de fino a grueso mientras que el algoritmo 
de Decodificación va de grueso a fino; ambos son algoritmos cuyo costo computacional es 
de O(Nq) operaciones ((Nq — iV¿) • (s + 1) sumas y (iVo — Nl) ■ s multiplicaciones). 

Notar además que los algoritmos de Codificación y Decodificación representan una 
transformada de ondelette exacta, pues u = M.^ 1 (Mu). 



2.2. Análisis de multiresolución por medias en celda 

En esta sección se considera la sucesión de iVo valores 

ü° = {«5}fi x (2.19) 

que se interpretarán como medias en celda (cell-averages) de cierta función u(x) sobre la 
malla fina G°: 

ü °j = T- í 3 u(x)dx, 1 < j < N . (2.20) 

ooooooooo 
k-1 I 1 1 1 1 1 1 1 1 





k | 1 1 1 

12 3 4 

Figura 2.2: Diferentes escalas de medias en celda 



Se consideran las mallas anidadas G h , 1 ^ k ^ L y se definen 



ü) = ^r \ 3 u(x)dx, 1 < j < N k . (2.21) 
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Se sigue de esta definición y de (2.5) que 



1 



hk Jx^_ l 
1 / 



u(x)dx 



'2j-l 



fc-i 

2¿-2 



■2j 



u{x)dx + / u{x)dx 



2j-l 



Por lo tanto 1 ^ j ^ iV fcj puede ser calculado en forma directa del dato inicial 

ü°, y sin ningún conocimiento explícito de la función u(x), mediante el algoritmo 



FOR k = l,2,...,L 
FOR j = l,...,N k 



(2.22) 



u 



fc-l^ 

2j 



Considerar la primitiva de u(x) 



U(x) = ¡ u(y)dy, 



(2.23) 



y observar que conocer las medias en celda u k es equivalente al conocimiento de los valores 
puntuales U k de la función primitiva, es decir, 

U k = {U k }f^ <-» ü k - s,-M^ 

lo cual es evidente de las siguientes dos relaciones: 

j 



U k = U{x k ) = P u{y)dy = V P u{y)dy = h k ¿ 
^ ¿=l^ti ¿=i 

_ fc _ C/(4) - U(aj_j) _ U k - U k _, 



u, 



hk hk 

~i nn T-inmVilí-i nnlnnlnv T 



(2.24) 



(2.25) 



En consecuencia conociendo los valores de u es posible calcular U y utilizar una función 
de interpolación para aproximar el valor ausente ^J-l' 1 ^ J ^ -^c por Ü k J~\, es decir, 



u 2j"-l ~~ - L V J '2j-V u 



(2.26) 



Con esto, y teniendo en cuenta que U k = U 2 - , es posible lograr una aproximación ü k 
para ü k ~ l mediante 



JT k ~ l - T~T k 

k-l _ U 2j-1 U j-1 -k-1 
' U 2j 



U 2j-l 



ifc-1 



ffk _ frk-1 
U j U 2j-1 

hk-l 



(2.27) 
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Notar que 



por lo tanto ü 2 ■ 1 puede calcularse a partir de u k y ü^j-i mediante 



ü*7 1 = 2üj-i^-_i 1 . (2.29) 



Se denota por d k (ü°) = {d k }f^ a la sucesión de errores de aproximación cometidos al 
fc-l \N 



predecir {tí^í-i}^! desde u k 



= vt^- 1 ^- 1 ^'^- 1 . (2.31) 



fc-l 

Análogamente al caso de valores puntuales, puede concluirse que existe una transformación 
uno a uno entre vP y su representación de multiresolución 

Ü M = (d 1 ,...,d L ,ü L f, (2.32) 

que se denota por 

ü M = Mü°, u° = M- l ü M . (2.33) 

En (2.31) el valor en x 2 j\ se ca l CUia a partir de la función polinomial que interpola 
los puntos (Uj_ 8 , . . . , U k +S _ 1 ). Utilizando lo visto anteriormente para el caso de valores 
puntuales, y como 2hk-i = h^, se obtiene 

2/ifc 



4 =fi -_?tí«(fcte^i, (2 . 34) 



con los Pi calculados en (2.16). 

De este modo, los coeficientes de ondelette están dados por 



-l 



= - ^ - E ^ ü "+i - ü ¿-i)> 1 < j < N k- ( 2 - 35 ) 

1=1 

Notar que se utiliza el esténcil (u k _ s+1 , . . . ,ü k +s _ l ) y por lo tanto el orden de precisión 
correspondiente es f = 2s — 1, y los coeficientes correspondientes 7/ son 

r = 3^ 71 = -1/8 
r = 5 => 71 = -22/128, 72 = 3/128 

Ver detalles en el apéndice A. 2. 
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Cuando se utiliza interpolación central (o cualquier interpolación independiente de los 
datos), se tiene que M es un operador lineal que puede ser expresado por una matriz de 
Nq x No. En el caso de que T(-,U k ) sea el especificado en la sección anterior, con r y 
s dados, las transformaciones en (2.32) pueden ser llevadas a cabo sin escribir la forma 
explícita de la matriz, calculadas mediante los algoritmos siguientes: 

Algoritmo de Codificación 



ü M = Mü L 



FOR k = l,2,...,L 



U; 



d) 



l 2j 



i k ~ l 



+/ 



1 < 3 < N k 



(2.36) 



Algoritmo de Decodificación 



i° = M~ l u M 



FOR k = L,L-1,...,1 
FOR j = l,...,N k 



j+i "./ / ' ' (l .r 



l 2j-l 



ü k + A, ü% 



A. 



(2.37) 



Ambos algoritmos poseen un costo computacional de O (No) operaciones ((Nq — iV¿) • 
(s + 2) sumas en ambos algoritmos y (A^o — N£) • s multiplicaciones en ambos algoritmos). 

Es interesante observar que dado que u° es equivalente a U°, también % es equivalente 
a Um, la representación de multiresolución de los valores puntuales de la función primitiva 

U(x) 

(d\d 2 , ...,d L , ü L ) T = ü M <-» U M = (D\D 2 ,..., D l , U l ) t . 
Además la transformación entre d k (vP) y D k (U°) está dada por 

d k (ü°) = D>;(U )/h k _ 1 . (2.38) 



2.3. Análisis de regularidad 

El análisis de multiresolución será de gran utilidad para obtener un algoritmo de com- 
presión de datos de las medias en celda. Luego se estudiará su aplicación a la solución 
numérica v n del esquema conservativo 

u n+l = u n_ A( j._j._ i)j X = r/h. (2.39) 
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Utilizando resultados de interpolación estándar y notando que U(x) es más suave que 
u(x), se obtiene de (2.38) la siguiente caracterización cualitativa de d k (ü°) (ver [26]): 



Teorema 1 Si la función u{x) en x = x* posee p — 1 derivadas continuas y una dis- 
continuidad de salto en la derivada p—ésima, entonces en los puntos x k cercanos a x* se 



tiene 



d k {ü 



~ < 



(h k )P[u^}, siO^p^f, 
{hk) p u^ p \ si p > r, 



(2.40) 



donde r es el orden de precisión de la aproximación (r = r — 1), 1 ¡/ [] denota el salto 
en la discontinuidad. 



Dem: Sea I(x, U ) como en (2.14). Se tiene que 



j+s-l 



U(x)=l(x,U k - l ) + U[x k r 1 s ,...,x l ;-l 1 ,x} ][ ( 



X — X A 



(2.41) 



%=]- 



con x S [x k _\,x k l ]. Notar que si u(x) tiene p — 1 derivadas continuas en x* y una 
discontinuidad de salto en cerca de x* , entonces U(x) tiene p derivadas continuas en 
x* y una discontinuidad de salto en U^ p+1 ^ cerca de x*. Con esto, de [1] se deduce que 



U[x\ 



,-fc-ii 
• • • ' x i+t 1 



0([1/(p+ 1 )]) . n , . . , 

k 

0(\\U®\\), sit<p+l. 



(2.42) 



Dado que D| = -X^*^, [/ fe ), l a relación (2.41) conduce a 



„k-l 



D k {U) = U[x k _ 



... ,x 



j+s-l 

_fc-l i TT /.Jk-l 



j+s-l > x 2j-lJ 



n 



v x 2j-l X ¿J' 



(2.43) 



y teniendo en cuenta que a?2j— i — x i es aproximadamente del orden de hk, con i £ {j 
s, . . . , j + s — 1}, se obtiene que 



D k AU) ~ <( 



(2.44) 



||C/ (r) ||/i£, sit<p + l. 
Finalmente, de (2.38), (2.44) y remarcando que (x) = w n >(x), se obtiene (2.40). 



□ 
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Ahora, la ecuación (2.40) en el nivel k — 1 corresponde a 



d 3 ~ < 



(/i fc _i)f[«W], «0<p<r-l, 
(^fc-i) r ~ 1 ^ r ~ 1 ^ sip>r—l, 



(2.45) 



y como /ifc = 2/ifc_i, entonces 



, 2"f(/i fc _ 2 )P[«(p)], «0<p<r-l, 
4 ~ i (2-46) 
2" r+1 ( / ifc-2) r_1 n (r - 1) , «p>r-l, 



Por lo tanto 

«2 _p |dJ|, p = mín(p,f). (2.47) 
Pueden obtenerse entonces algunas conclusiones útiles 

■ Lejos de las discontinuidades, los coeficientes dj decrecen a medida que se va a niveles 
más finos. 

■ La tasa de decaimiento de los coeficientes dj es determinada por la regularidad local 
de la función y el orden de precisión de la aproximación. 

■ En la vecindad de una irregularidad de u(x), los coeficientes dj permanecen del 
mismo orden 0([i¿]), independiente del nivel de refinamiento. 

Por lo tanto el análisis de multiresolución de ü puede verse como un estudio de la regu- 
laridad local de u(x). 

Puede hacerse un análisis de regularidad similar si se considera el caso de valores pun- 
tuales en vez de medias en celdas. De forma análoga, Kozakevicius (ver [29]) propone que 
dependiendo de la regularidad de la función, un gran número de coeficientes de ondelet- 
te pueden ser extremadamente pequeños, y por lo tanto podrían ser descartados de la 
representación de multiresolución. 



2.4. Compresión de datos 

La idea principal es reducir la cantidad de datos mediante una técnica de truncamiento, 
que consiste en hacer ceros los coeficientes que están por debajo de una tolerancia prescrita. 
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Sea tr £k el operador de truncamiento definido por 



4 = tr £k (4) = ^ 



v 



0, si |dj| < e k 
dj, en otro caso. 



(2.48) 



Sea u"m el resultado de la operación de truncamiento aplicada a um 



um = {d l ,d 2 , . . . ,d L ,u L ). 



(2.49) 



Si se aplica el algoritmo de decodificación al dato truncado u"m, se obtiene una aproxima- 
ción u° = M _1 um, que por [26] se sabe que permanece cerca del dato inicial u°. 

Dado que se está en el caso de multiresolución por valores puntuales de u, cada coefi- 
ciente de ondelette está relacionado con una posición específica en la malla fina uniforme 
y por lo tanto los procesos de codificación y decodificación pueden ser simplificados. Los 
coeficientes dj se calculan entonces sólo para decidir si x k - seguirá o no en la malla y se evita 
así construir la representación de multiresolución completa [29]. Esto quiere decir, que en 
estos puntos, la información sobre la función puede ser obtenida mediante interpolación. 

La representación de u° al cabo de este proceso, contendrá sólo los valores puntuales en 
las posiciones asociadas a coeficientes de ondelette significativos, y los puntos en el nivel 
más grueso. Esto se conoce como representación puntual esparsa de u, y se denota por us- 

La elección de e k puede variar de acuerdo a las propiedades de los espacios funcionales 
¡ ], o suavidad de la función [26]. En este caso, con e fijo, los niveles de tolerancia en 
cada nivel estarán dados por = e/2 L ~ k . Notar que a escalas más finas, es más 
pequeño; esto con el fin de preservar la información asociada a la parte regular del dato 
inicial y descartar perturbaciones de alta frecuencia (pues una señal regular posee mayores 
coeficientes de ondelette en escalas más gruesas y una señal perturbada, o una función con 
singularidades, posee mayores coeficientes de ondelette en escalas más finas). Además esta 
elección de es óptima en el sentido que mantiene la mejor relación entre compresión de 
datos y disipación de información durante la evolución temporal de la solución. 

La representación puntual esparsa u$ también incluirá algunos safety points necesarios 
para evitar la disipación numérica; este corresponde al operador de extensión E. Los 
safety points serán incluidos en las vecindades de puntos cuyos coeficientes de ondelette 
son significativos [29]. Se incluirán dos tipos de safety points: Puntos en el mismo nivel 
de multiresolución que el coeficiente de ondelette respectivo (con el fin de mantener la 
calidad del transporte de información desde un punto a su vecino en la malla) y puntos 
en un nivel de multiresolución más fino que el nivel del coeficiente de ondelette (sólo si el 
detalle es mayor que una tolerancia adicional 2e¿., esto con el fin de mejorar la captura de 
choques). 
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Figura 2.3: Secuencia de operaciones para obtener la representación puntual esparsa de una fun- 
ción. DWT: transformada del dato inicial, tr £fc : operador de truncamiento, E: inclusión de safety 
points, IWT: transformada de ondelette inversa y R: reconstrucción de malla uniforme. 



2.5. Estructura de datos 



Dado las características de los problemas hiperbólicos que poseen discontinuidades que 
se propagan, el número de puntos en la representación puntual esparsa es mucho menor 
que el número de puntos en la malla fina uniforme. Luego, será de gran utilidad almacenar 
la información relevante en algún tipo de estructura que saque provecho de ello, tal como 
se hace en [29] (MORSE, SPARSE, etc.) 



1 2 3 


4 5 6 7 .... 


21 22 ... 128 


• • 


mm i m 












Vector de valores puntuales 


r ir ir ir ir 

L 








1 2 3 4 5 



Vector de posiciones 2 | 6 1 2l| 22j i2s| 



1 2 3 4 5 
J Valores puntuales en posiciones significativas 

~] Posiciones de coeficientes descartados 

Figura 2.4: Ejemplo de almacenamiento de datos sólo para posiciones significativas de la repre- 
sentación truncada (MORSE o SPARSE). 
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Capítulo 3 

Caso hiperbólico 



En esta sección se presenta una forma eficiente de resolver leyes de conservación hi- 
perbólicas mediante un método de alta resolución en mallas generadas por ondelettes 
desarrollado por Kozakevicius y Santos [29]. La eficiencia de este método se basa en la 
asociación de dos técnicas independientes: mallas adaptativas generadas por una transfor- 
mación de ondelettes [26, 14, 28] y métodos de alta resolución basados en interpolaciones 
ENO para el cálculo de los flujos [33, 29]. 



3.1. Esquema ENO Lax-Friedrichs 

Se necesitan esquemas conservativos para la parte espacial del operador (forma semi- 
discreta) 

^KW) = ^t(w-^i/ 2 ), 

donde /j+1/2 = f(uj- r ,---,Uj- s ) es el flujo numérico, en que la primera posición del 
esténcil j — r es elegida mediante un algoritmo ENO, manteniendo la relación j — r < 
j + 1/2 < j — s. Esta función de flujo numérico es Lipschitz continua en sus argumentos 
y es consistente con el flujo exacto, es decir, f(u, . . . , u) = f(u). 

Para lograr un alto orden de aproximación para j¿ , se utilizarán posiciones escalonadas 
auxiliares {xj +1 / 2 }j [ ] con respecto a la malla gruesa esparsa. El flujo numérico evaluado 
en estas posiciones se obtiene mediante interpolación ENO. 

Es necesario considerar esquemas upwind en la construcción del flujo numérico con 
el objetivo de mantener la estabilidad del esquema. Para ello se utilizará la forma más 
sencilla, más robusta y menos costosa de obtener esquemas upwind sin violar condiciones 
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h. 

i 




Figura 3.1: Componentes de la separación del flujo numérico en la frontera, hj es el inter- 
polador ENO para la celda ]xj_ l / 2 , Xj + i/ 2 [ Y hj+i es e l interpolador ENO para la celda 
]x j+1 / 2 ,x j+1+1 / 2 [. 



de entropía de la solución. Esta es, la separación de flujo de Lax-Friedrichs: 

f(u) = f + (u) + f-(u), f+( u )= l -{f{u) + au), r{u) = \{f{u)-au), 

donde 

a = máx|/'(u)|. (3.1) 

u 

El número de puntos escogidos para la reconstrucción depende del orden de la interpola- 
ción. En este caso, se utilizará interpolación cúbica. 

El flujo numérico en las posiciones de la malla auxiliar corresponde a la suma de las 
aproximaciones generadas para cada parte de la separación de flujos 

fj+i/2 = í j +1/2 + fj+i/r ( 3 - 2 ) 

Notar que fj~ + i/ 2 Y f j+1/2 son a P rox i m aciones para el mismo borde Xj + i/ 2 del volumen de 
control ]xj_x/2,Xj + i/2[, obtenidas mediante interpoladores distintos. 

Notar además, que una vez que se elige el número de puntos en el esténcil, este per- 
manece igual para todos los puntos de la malla G k . Esta forma de construir predictores 
para la transformada de ondelette no considera la suavidad local de la función a ser inter- 
polada. Si la función es suave a trozos, una aproximación de esténcil fijo puede compor- 
tarse inadecuadamente cerca de las irregularidades, generando oscilaciones en las celdas 
correspondientes. Estas oscilaciones (conocidas como Fenómeno de Gibbs en métodos es- 
pectrales) ocurren debido a que los esténciles contienen una celda discontinua (volumen 
de control que contiene una irregularidad), es decir, poseen un punto Xj bastante cerca de 
una irregularidad. Además, cada vez que el esténcil cruza una singularidad, la calidad de 
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la interpolación se ve reducida. Cuanto mayor es el grado del interpolador, mayor es la 
región afectada por la singularidad. 

La idea es entonces utilizar interpolación ENO {Essentially non oscillatory), que au- 
menta la región de precisión para el interpolador, eligiendo un esténcil diferente, para 
evitar las oscilaciones cerca de las discontinuidades. 

Se presenta a continuación la forma en que se prepara la reconstrucción ENO. Inicial- 
mente se conocen los valores de los flujos en la malla esparsa S. Se define V(xj + i/ 2 ), la 
primitiva de la componente de separación de flujo en la malla auxiliar con respecto a S. 
Se construirá un polinomio interpolador por partes de V, en la variable x: H(x, V), sobre 
la malla auxiliar, es decir, 

i 

H(x j+1/2 , V) = V j+ i/2 = V{x j+l/2 ) = f(x h ), 

H(x, V) = q m (x, V), Xj-1/2 < x < x j+í / 2 , 

donde q m es el único polinomio interpolador de grado m, que utiliza m + 1 puntos conse- 
cutivos (x im(i) , . . . ,x im(i)+m ), incluyendo a Xj_ 1/2 y x j+1/2 . 

Notar que dependiendo de la elección del primer punto del esténcil i m (j), existen m 
polinomios interpoladores posibles. ¿Cuál elegir? El esténcil asociado a [xj-i/ 2 ,Xj +í / 2 ] 
será aquel tal que V{x) es más suave (en un sentido asintótico) y el valor x donde se 
evaluará el interpolador, será Xj_y 2 o Xj +1 / 2 . 

La información de la suavidad de V puede obtenerse de las diferencias divididas: 

= V(xj_i/ 2 ) 

, , w[x j -i/2 + l,...,X j -U2 +k ]-w[x j -U2,...,Xj-i/2 +k - 1 ] 

w[x J _ l/2 ,...,x j _ 1/2+k ] = . 

x j-l/ 2 +k — x j-l/2 

El siguiente Teorema (ver [29, 17]) entrega un criterio para medir asintóticamente la 
suavidad de una función. 

Teorema 2 Si f(x) es C°°([xí, Sj+jj.]), entonces 

1 d k 

w[xí, . . . ,x i+k ] = — -j-¡¿f{€i,k), Xi < £ i¡k < x i+k ; (3.3) 
pero si f(x) tiene una discontinuidad de salto en sup—ésima derivada, ^ p ^ k, entonces 
w[x i ,...,x l+k ] = 0(d- k+ n[f {p) ], d=\x l+k - Xi \. (3.4) 

Luego, utilizando \w[xj_i/ 2 , . . . , Xj_i/ 2+k ]\ es posible medir asintóticamente la suavidad 
de f(x) en [xj_i/ 2 , Xj_i/ 2+k \: El mejor esténcil, será aquel asociado a la diferencia dividida 
más pequeña [17]. 
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La cuestión es ahora, cómo hallar i m (j). Para ello se seguirá el siguiente procedimiento 
([29]): 

1. = j — 2j donde q\ es el polinomio interpolador para V en £j_i/2 y x j+i/2- 

2. Suponer un polinomio interpolador de grado ra, q n para V en • • • , £j n (j) +n . 

3. De acuerdo con la diferencia dividida más pequeña, q n +\ comenzará con i n+ i{j) = 
in(j) — 1 (si el siguiente punto elegido está a la izquierda del último punto en el 
esténcil) o con í n+ i(j) = i n (j) (si el siguiente punto elegido está a la derecha del 
último punto en el esténcil). 

Las aproximaciones hj_i/-¿ y para cada componente de la separación de flujo serán 

entonces la derivada de H evaluada en Xj_u 2 y Xj+1/2 respectivamente. Una vez calculado 
el flujo numérico en las posiciones auxiliares, se obtiene una aproximación de alto orden 
para el término de la derivada espacial en las posiciones de la malla esparsa S. 

3.2. Evolución temporal 

Notar que se está frente a un proceso de discretización en dos etapas, primero se ha 
discretizado sólo espacio, dejando el problema continuo en tiempo. Esto conduce a las 
llamadas ecuaciones semi- discretas. La discretización puede hacerse utilizando un método 
numérico estándar para sistemas de ecuaciones diferenciales ordinarias. Este mecanismo 
es particularmente ventajoso en el desarrollo de métodos con orden de precisión mayor a 
dos, ya que permite alcanzar de forma relativamente sencilla la misma precisión espacial 
y temporal. 

Los experimentos realizados en [33] indican que las formulaciones semi-discretas con 
discretización temporal Runge-Kutta TVD desarrollados por Shu y Osher no generan os- 
cilaciones para CFL ^ 0,5 aproximadamente, y son óptimas en el sentido de que permiten 
el mayor CFL para esquemas explícitos, CFL = 1. 

Se utilizarán entonces métodos Runge-Kutta TVD de segundo o tercer orden. 

R-K TVD óptimo de segundo orden: 

= u n + AtC{u n ) 
u n+1 = ] ¿ u n + ^u (1) +At£(u w ), 
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R-K TVD óptimo de tercer orden orden: 

u (1) = u n + AtC{u n ) 

= -u n + \u^ + AtC{u^) 
4 4 

u n+1 = ~u n + 2 uW + 2 AtC(uW), 
3 3 3 V ' 

con C{u) = -{Ax)-\fj{u) - f á -i(v)). 

La alta resolución (asociada a discretizaciones espaciales TVB, ENO o TVD) es necesa- 
ria para asegurar estabilidad. En los pasos intermedios del esquema de evolución temporal, 
se conserva la malla esparsa del paso n. 

3.3. Adaptatividad de la representación esparsa 

Con el fin de actualizar la malla esparsa, es necesario aplicar el operador de reconstruc- 
ción R para reconstruir la solución en la malla uniforme. Una vez aplicada la transformada 
de ondelette, el operador de truncamiento y el operador de extensión, puede llevarse a cabo 
la evolución temporal. 

Dado que recalcular la malla es costoso, puede utilizarse la misma malla para varios 
pasos temporales. Para problemas donde la velocidad de la onda es baja (en el sentido 
CFL), es posible utilizar la misma representación puntual esparsa para 5 o más pasos 
temporales sin aumentar la disipación numérica, y luego realizar la actualización de la 
configuración. Para problemas con una alta velocidad de onda, la reconstrucción de la 
malla puede hacerse cada dos pasos temporales, sin afectar la calidad de la solución [29]. 

Cuando se trabaja con ecuaciones multivariadas se construye una malla esparsa "uni- 
ficada". Es la unión de las posiciones significativas de la representación esparsa de cada 
componente y todos los safety points necesarios para la evolución, en cada componente. 
El criterio para la malla unificada es bastante simple. Una vez que una posición tiene 
asociado un coeficiente de ondelette significativo en cualquier componente del vector de 
cantidades, tal posición debe permanecer en la malla unificada, y todas las componentes 
del vector de cantidades deben tener sus valores puntuales en esta posición. Lo mismo 
sucede con el operador de extensión. 

Notar que como cada variable del vector de cantidades desarrolla discontinuidades 
bastante localizadas, la malla unificada seguirá siendo esparsa [26, 27, 29]. 

La actualización de la malla es análoga al caso escalar. Los mismos operadores de- 
ben ser aplicados a cada componente del vector de cantidades para obtener la siguiente 
configuración de la malla unificada y realizar la evolución temporal. 
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3.4. Método adaptativo de alta resolución 

Dado el número de puntos en la malla fina, Nq, el número de niveles de multiresolución, 
L, el grado r del predictor intermallas y del interpolador ENO, el nivel de truncamiento 
dadas además las condiciones de contorno e inicial de la ley de conservación, el algoritmo 
del método descrito puede ser resumido como sigue: 

1. ■ Transformada de ondelette discreta {DWT) (u operador de codificación M) 

aplicada al dato inicial. 

■ Representación puntual esparsa {SPR) de la solución. Esta incluye trunca- 
miento, extensión, y transformada inversa de ondelette (IWT) (u operador de 
decodificación M _1 ). 

2. ■ Cálculo del flujo exacto en malla esparsa (correspondiente al nivel más fino de 

multiresolución) . 

■ Cálculo del valor global de a (3.1). 

■ Cálculo de Ai para la evolución temporal: Ai = CF ^' h o ; donde ho es el paso 
espacial en la malla fina. 

■ Factorización Lax-Friedrichs del flujo exacto: f + yf~- 

■ Cálculo del flujo numérico fj+i/2- 

• para / + , construir la aproximación ENO h~ . 

• para /", construir la aproximación ENO h + . 

• f j+l/2 = f 3 + +1 / 2 + fj+1/2' 

3. ■ Evolución temporal: Runge-Kutta TVD de segundo o tercer orden. Se necesitan 

pasos intermedios. 

■ Repetir 2. para la solución intermedia necesaria para 3. 

■ Evolución temporal de la solución intermedia (El método Runge-Kutta TVD 
de segundo orden completa el paso temporal, el método Runge-Kutta TVD de 
tercer orden necesita otro paso intermedio). 

4. Aplicación del operador de reconstrucción de la solución en malla fina R. 

5. Volver al., aplicar DWT a la solución obtenida y repetir (ver [29]). 
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3.5. Resultados numéricos 



En esta sección se reproducirán algunos resultados obtenidos por Harten [26]. Para ello 
se aplicará el algoritmo de multiresolución a la solución numérica de una ley de conserva- 
ción, tomando como modelo la ecuación de Burgers (caso escalar y unidimensional) 

u t + {u 2 /2) x = (3.5) 

asociada a la condición inicial 

1, si \x\ < 1/2 



u(x, 0) 



(3.6) 

0, si 1/2 < \x\ < 1. 



Se utilizan condiciones periódicas en i = - 1 y i = 1. Se opera hasta antes de que las 
discontinuidades alcancen las fronteras del dominio. 

El primer objetivo es mostrar la relación existente entre la capacidad de compresión de 
este método de multiresolución y las propiedades de aproximación de las técnicas de re- 
construcción utilizadas. La localización de los coeficientes de ondelette que están por sobre 
una tolerancia prescrita, ayuda a visualizar esta conexión. Recordar que los coeficientes 
de ondelette dj representan los errores cometidos en el proceso de predicción y están di- 
rectamente relacionados a errores de interpolación, los cuales son pequeños en regiones de 
suavidad. En las proximidades de las singularidades el proceso de reconstrucción podría 
conducir a regiones de exactitud pobre, por lo tanto, se examina el efecto del esquema 
de compresión basado en la multiresolución. Como una medida de la mejora en velocidad 
alcanzada mediante la utilización del análisis de multiresolución, se presenta la tasa de 
compresión o eficiencia \x [3, 26] definida por ii = jVo/P+ji3"T ' donde D n es el conjunto de 
coeficientes de ondelette significativos, en todos los niveles de multiresolución, en el paso 
temporal ra. 

Las figuras 3.2 a 3.9 y las tablas 3.1 y 3.2 resumen el resultado de los test numéricos 
realizados. En cada figura, la parte izquierda representa a la solución numérica con asteris- 
cos. La parte derecha muestra el conjunto de los coeficientes de ondelette significativos en 
el plano x — k, dibujando un + alrededor de cada (xj, k). Cada tabla muestra resultados 
de multiresolución para la solución numérica de la ecuación de Burgers para diferentes 
tiempos t. Se muestra la tasa de compresión ri, proporción V (entre el tiempo total de 
CPU de la solución numérica sin multiresolución y el tiempo total de CPU de la solución 
numérica con multiresolución) y los errores e p = \\u n — u^ IR \\ p , p = 1,2, oo, donde 

eoo = máx \u] - umrj \, 1 < 3 < 

y 

i/p 
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En ambas tablas se verá que el error obtenido es menor que la tolerancia prescrita. La 
norma £00 obtiene el menor error principalmente en funciones discontinuas [26]. Es im- 
portante precisar que los errores son calculados entre la representación puntual esparsa 
y la solución en malla fina, aún cuando la longitud de estos vectores no coincide (ya que 
existen posiciones en la malla fina para los cuales no corresponde ningún punto en la 
representación puntual esparsa). 

Se presentan los resultados correspondientes para el caso de No = 257 puntos en la 
malla fina con L = 7 niveles de multiresolución y el caso de iVo = 1025 puntos en la 
malla fina con L = 10 niveles de multiresolución. En ambos casos se utiliza una tolerancia 
de truncamiento = e/2 L ~ k , condición CFL = 0,5, multiresolución con interpolador 
cúadrático, flujos numéricos calculados mediante reconstrucción ENO de segundo orden 
(ver sección 4.1.5) y evolución temporal Runge-Kutta de orden 2 (4.7). 



t 


V 




ei 


e-2 








0.16 


1.9330 


19.7633 


8.89xl0~ 7 


1.92x10" 


-5 


1.80x10- 


-4 


0.47 


1.8334 


19.8122 


1.99xl0~ 6 


3.15x10" 


-5 


6.14x10" 


-5 


0.62 


1.7696 


19.4591 


2.46xl0~ 5 


3.58x10- 


-5 


5.91x10" 


-5 


0.78 


1.6881 


19.7633 


2.92xl0~ 5 


3.96x10- 


-5 


5.77x10- 


-5 



Cuadro 3.1: Solución numérica de la Ecuación de Burgers, condición inicial (3.6). Toleran- 
cia prescrita e = 10~ 5 , Nq = 257 puntos en la malla fina y L = 7 niveles de multiresolución. 



t 


V 




ei 


e2 








0.16 


2.7872 


53.9446 


8.94xl0~ 6 


1.99x10- 


-6 


4.79x10- 


-6 


0.47 


2.5986 


53.0172 


2.09xl0~ 5 


3.01x10- 


-6 


5.59x10- 


-6 


0.62 


2.6170 


53.5019 


2.49xl0~ 5 


3.99x10- 


-6 


7.26x10- 


-6 


0.78 


2.5029 


53.2874 


2.97xl0~ 5 


4.26x10- 


-6 


1.88x10- 


-5 



Cuadro 3.2: Solución numérica de la Ecuación de Burgers, condición inicial (3.6). To- 
lerancia prescrita e = 10~ 3 , iVo = 1025 puntos en la malla fina y L = 10 niveles de 
multiresolución . 
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Ecuación de Burgers en 1D, CFL=0.4, t=0.1 6 




-0.5 



0.5 



Coeficientes de ondeletl significativos 

7 1 1 1 



o 5 



S> 3 



-0.5 0.5 

x 



Figura 3.2: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,16 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con e — 10~ 5 , 
A<o = 257 y L — 7. Derecha: Estructura de coeficientes de ondelette significativos correspondientes. 



Ecuación de Burgers en 1D, CFI_=0.4, t=0.47 




Coeficientes de ondelett significativos 

7 1 1 1 



o 5- 
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2 - 
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-0.5 
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Figura 3.3: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,47 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con 
e = 10~ D , A^q = 257 y L = 7. Derecha: Estructura de coeficientes de ondelette significativos. 
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Ecuación de Burgers en 1 D, CFL=0.4, t=0.62 Coeficientes de ondelett significativos 

1 1 1 1 7 1 1 1 




Figura 3.4: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,62 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con 
e = 10~ a , No = 257 y L — 7. Derecha: Estructura de coeficientes de ondelette significativos. 



Ecuación de Burgers en 1 D, CFL=0.4, t=0.78 Coeficientes de ondelett significativos 

1 1 1 1 7 1 . 1 




i , , , 1 1 1 , , , 1 

-1 -0.5 0.5 1 -1 -0.5 0.5 1 

x x 

Figura 3.5: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,78 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con 
e = 10~ D , N = 257 y L — 7. Derecha: Estructura de coeficientes de ondelette significativos. 
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Ecuación de Burgers en 1D, CFI_=0.4, t=0.1 6 



Coeficientes de ondelett significativos 
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Figura 3.6: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,16 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con 
A<o = 1025, L = 10, e = 10~ 3 . Derecha: Estructura de coeficientes de ondelette significativos. 



Ecuación de Burgers en 1D, CFL=0.4, t=0.47 
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Coeficientes de ondelett significativos 
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Figura 3.7: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,47 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con 
A^o = 1025, L = 10, e = 10~ 3 . Derecha: Estructura de coeficientes de ondelette significativos. 
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Ecuación de Burgers en 1 D, CFL=0.4, t=0.62 Coeficientes de ondelett significativos 

1 1 1 1 10 i 1 1 1 




I 1 1 1 1 I 1 1 1 1 

-1 -0.5 0.5 1 -1 -0.5 0.5 1 

x x 



Figura 3.8: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,62 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con 
A<o = 1025, L — 10, e = 10~ 3 . Derecha: Estructura de coeficientes de ondelette significativos. 



Ecuación de Burgers en 1 D, CFL=0.4, t=0.78 Coeficientes de ondelett significativos 

1 , , 1 10 | , 1 , 




Figura 3.9: Izquierda: Solución inicial (rayas) y solución numérica de multiresolución (asteriscos) 
en el tiempo t = 0,78 para la ec. de Burgers en ID asociada a la condición inicial (3.6), con 
No = 1025, L = 10, e = 10~ 3 . Derecha: Estructura de coeficientes de ondelette significativos. 
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Capítulo 4 

Caso parabólico 



En este capítulo se aplicará el algoritmo de multiresolución a ecuaciones parabólicas. 
Se reproducirán los experimentos numéricos realizados por Roussel et al. [32] , Bihari [3] , 
Liandrat y Tchamitchian[30]. 

4.1. Método numérico 

A continuación se presenta un método general de volúmenes finitos para ecuaciones 
hiperbólicas, incluyendo la descripción de los esquemas utilizados para la discretización 
espacial y evolución temporal [32]. 

4.1.1. Leyes de conservación parabólicas 

Se considera el problema de valores iniciales para una ecuación parabólica en (ce, t) £ 
ü x [0, oo[, ü C M d de la forma 



asociada a condiciones de borde apropiadas. 

Se considerará la restricción al caso en que el flujo difusivo se define por un operador 
gradiente, suponiendo difusividad constante v > 0, es decir, 




(4.1) 



u(x, 0) = uq(x) 



F(u,Vu) = f(u) - vVu. 
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Para la ecuación de convección-difusión en ID, se tiene (c > 0) 

f{u) = cu, 
S(u) = 0. 

En el caso de la ecuación viscosa de Burgers en ID, se tiene 

u 2 

/(«) = y- 

S(u) = 0, 

y para la ecuación de reacción-difusión (a > 0,/3 > 0), 

/(«) = o, 

% = v( 1-u ) exp ^i — s — r 

2 a(l — uj — 1 

Se define el término fuente y de divergencia por 

T>(u, Vu) = -V • F(u, Vm) + S(u). 
Luego (4.1) puede escribirse como 

du 



()l V{u,Vu). (4.2) 



4.1.2. Discretización 



Para discretizar (4.2), se utiliza una formulación de volúmenes finitos en la forma 
conservativa estándar. En el caso general, considérese el dominio computacional íl y una 
partición de él en volúmenes de control (0¿)¿ 6 Aj A = {1, . . . ,i ma x}- Se denota entonces 
por (¡i{t) al promedio de cierta cantidad q sobre íí¿ en el instante t, 

?<(*) = 7777 / ?0M)<&- (4-3) 



Integrando (4.2) y promediando sobre ííj, 



1 í du 1 

(x,t)dx = -—l V(u(x,t),Vu(x,t))dx. (4.4) 



Luego 

^(t)=A(t). (4-5) 
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Si se aplica el teorema de la divergencia, se obtiene 

Vi = —^-;¡ F(u,Vu)-ai(x)dx + Si(t), (4.6) 

donde <r¿(x) es el vector normal hacia íí¿. La conservatividad en el cálculo del flujo se 
garantiza si y sólo si, para dos volúmenes de control adyacentes íí^ y íí¿ 2 , el flujo que va 
de ííjj a Q¿ 2 se equilibra con el flujo que va de Í2¿ 2 a fi^ . 



4.1.3. Integración temporal 

Notar que se está frente a un proceso de discretización en dos etapas, debido a la 
adaptatividad de la discretización espacial. Primero se discretiza sólo en tiempo, y luego 
en espacio. Esto conduce a las ya mencionadas ecuaciones semi-discretas (ver sección 3.2). 
La discretización puede hacerse utilizando un método numérico estándar para sistemas 
de ecuaciones diferenciales ordinarias. Este mecanismo es particularmente ventajoso en 
el desarrollo de métodos con orden de precisión mayor a dos, ya que permite alcanzar 
de forma relativamente sencilla la misma precisión espacial y temporal. Los experimentos 
realizados en [33] indican que las formulaciones semi-discretas con discretización temporal 
Runge-Kutta TVD desarrollados por Shu y Osher no generan oscilaciones para CFL ^ 
0,5 aproximadamente, y son óptimas en el sentido de que permiten la mayor CFL para 
esquemas explícitos, CFL = 1 [33, 13, 29]. 

Se utilizará entonces un método explícito Runge-Kutta TVD de segundo orden que en 
este caso se expresa por 



_n+l/2 = Ü n + Atf) n 



"i 

= ~ 



2 

Notar que (4.7) también se conoce como Método de Heun [3]. 



(4.7) 



Si se denota por u n al vector (i¿™)¿ g A, entonces el operador de evolución temporal 
discreto E(At) está definido por 

Ü n+l = É(At) ■ ü n , (4.8) 

donde 

É{At) = I + ^ [V + V(I + AW)] . (4.9) 
La discretización del operador T> se describe en la siguiente sección. 
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4.1.4. Flujo numérico 

Considérese ahora un tiempo fijo t n . Para el caso unidimensional general, íí¿ es el 
intervalo [x¿_i/2i ^¿+1/2] de longitud Ax¿ = x i+ i/ 2 — %í-i/2- Mediante una discretización 
de volúmenes finitos estándar, la ecuación (4.6) puede escribirse como 

V* = - 1 ?-(F í+1 -F t A+S l , (4.10) 



2 



donde 



!+ 2 

con Ax ¿+ i = t;(Ax¿ + Axj+x)- El término f R denota, para la parte advectiva, la solución 
aproximada de Roe para el problema de Riemann [21], dados los estados de derecha e 
izquierda de u. La versión escalar correspondiente es 

f R (u~,u+) = l -[f{u-) + f{u + ) - \a(u~,u + )\(u + - «-)], (4.12) 



donde 

a(u~ , u + ) 



ñ u+)-f(u-) úu +^ u - 
u + —u ' ' ' 



f'(u + ), SÍ U + = U . 

Los valores de izquierda y derecha u~ 1 y ; respectivamente, son obtenidos mediante 
interpolación ENO de segundo orden (ver sección 4.1.5). 

Notar de (4.11) que los términos advectivo y difusivo son aproximados de diferente 
forma. Para la parte advectiva, se utiliza el esquema de Roe clásico con una interpolación 
ENO de segundo orden; mientras que para la parte difusiva, se escoge un esquema centrado 
en Ui de segundo orden. 

En [3] se prueba que el esquema global resultante, que es no lineal, 

* ~ ¿ (<* fe^) - <- («M-*») - ^' A y- ) + *. («3) 

es de segundo orden (en espacio). 

El término fuente es aproximado por Si ~ S(üi). Para un término fuente no lineal, esta 
elección también implica una precisión de orden dos [32]. 



4.1.5. Reconstrucción ENO de segundo orden 

Para obtener los valores de la función u en las fronteras de los volúmenes de control, 
se utiliza una reconstrucción lineal a trozos de u a partir de los valores de las medias en 
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celda. Es decir, los términos de izquierda y derecha u x y uT + x , respectivamente, son 

obtenidos mediante interpolación ENO de segundo orden [29, 33, 13, 32]. Este tipo de 
métodos utiliza una construcción adaptativa del esténcil a fin de evitar la generación de 
oscilaciones espúreas cerca de las discontinuidades. Se puede generar oscilaciones, pero 
del orden del error local de truncamiento en la parte suave de la solución. En este caso 
particular, se tiene 

ü~x = üi + \m (ü i+1 - üí,üí - üj-i) , (4.14) 

¿_1 " 2 ¿ 

üf.í = ü i+ i + -M(ü i+2 - ü i+1 ,ü i+ i - Ui) , (4.15) 

¿ T 2 ¿ 

donde M es el limitador Min-Mod, que escoge la pendiente mínima entre los extremos 
izquierdo y derecho, es decir, 



M(a,b) = < 



a, si \a\ ^ \b\, 
6, si \a\ > \b\. 



Notar que (4.11) es la forma semi- discreta de (4.8). (4.11) se resuelve utilizando una 
actualización temporal Runge-Kutta de segundo orden; por lo tanto se obtiene un esquema 
de segundo orden tanto en tiempo como en espacio. 

Mediante un argumento de producto tensorial, puede llevarse a cabo la extensión na- 
tural de la reconstrucción a 2D y 3D en geometrías cartesianas [32]. 

4.1.6. Solución exacta de la onda viajera 

Para formar una idea cualitativa de la estructura del choque, considérese la solución 
u(x,t) = u(ip), ip = (x — st)/v del problema de la onda viajera 

u t + f{u) x = vu xx , (4.16) 
( 

ul, si x < 0, 



u(x, 0) 



(4.17) 



ur < ul, si x ^ 0. 
La ecuación diferencial ordinaria resultante en tp puede integrarse para obtener 

- su + f(u) + c = u, (4.18) 
donde s y c pueden ser determinadas de las "condiciones de borde" 

lím u(ip) = ul, lira u(ip) = ur 

ip^—oo ip— >oo 
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como sigue 



c = su L -f(u L ), (4.19) 
s = (4.20) 

donde la velocidad de la onda s puede ser identificada como la velocidad del choque (asume 
la misma expresión que en el caso puramente hiperbólico). Una nueva integración de (4.18) 
entrega una fórmula implícita para u: 

''" V + ci. (4.21) 



f(u) — su + c 

En el caso particular de la ecuación de Burgers viscosa, es decir, f(u) = ¿i¿ 2 , se obtiene 



u(i/>) = u L tanh UL — x¡). (4.22) 



Ver detalles en [3]. 



4.1.7. Estabilidad numérica 

Como el paso temporal es el mismo para todas las escalas de multiresolución, la condi- 
ción de estabilidad es la correspondiente al esquema de volúmenes finitos en la malla fina. 
Si denotamos por Ax al menor paso espacial, el número CFL o está dado por 

Ai 

o- = u máx — . (4.23) 
Ax 

Para el caso lineal (ecuación de convección-difusión), si c es la velocidad, 

y el número de Reynolds Re está dado por 

Re = — . (4.25) 
v 

En [3] y [21] se muestra que una condición suficiente para asegurar la estabilidad del 
esquema de volúmenes finitos es 

a ^ rain — , — - . (4.26) 



V 2 ' Re 

Aún más, una condición suficiente para que el esquema sea TVD (ver apéndice B), es 

Re 

a ^ -J^. (4.27) 
fíe + 4 v ' 
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La mayor ventaja de utilizar un esquema explícito para el término difusivo, es que no 
se necesita resolver un sistema lineal. Sin embargo, esto generalmente implica que At = 
0(Ax 2 ). Sólo para el caso Re » 1 se puede esperar Ai = 0(Ax) [20]. 

A continuación se analizará un Esquema de multiresolución conservativo completamente 
adaptativo diseñado por Roussel et al. [32]. 

4.1.8. Árbol graduado dinámico 

El principio del análisis de multiresolución es representar un conjunto de datos dados en 
malla fina como valores en la malla más gruesa y un conjunto de detalles a diferentes escalas 
de mallas anidadas. Se propone organizar la estructura de datos como un árbol graduado 
dinámico, que posee una capacidad mayor de compresión que la estructura MORSE o 
SPARSE de la representación puntual esparsa. 

En la terminología de las ondelettes, una estructura de árbol graduado corresponde a 
una aproximación adaptativa en la que está garantizada la conectividad para la estructura 
de árbol. 

Para definir la estructura de árbol, se introduce la terminología utilizada por Cohén 

[14, 32] : 

■ La raíz es la base del árbol. 

■ Un nodo es un elemento del árbol. Cada volumen de control será considerado un 
nodo. 

■ Un nodo padre tiene 2 nodos hijos; los nodos hijos de un mismo nodo padre son 
llamados hermanos. 

■ Un nodo tiene vecinos cercanos en cada dirección, llamados primos cercanos. Los 
nodos hermanos son también considerados como primos cercanos. 

■ Un nodo es llamado hoja cuando no tiene hijos. 

■ Para calcular los flujos entrantes y salientes de cada hoja, se necesitan los primos 
cercanos. Cuando alguno de ellos no existe, se crea una hoja virtual (representada 
por rayas en la figura 4.1). Esta no se considera como un nodo existente, sino sólo 
se utiliza para calcular flujos. 

Un árbol dinámico es un árbol que cambia en el tiempo. Si es necesario, algunos nodos 
pueden ser agregados o quitados. Para permanecer graduado, el árbol debe respetar las 
condiciones siguientes: 
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■ Cuando un hijo es creado, todos sus hermanos son creados en el mismo tiempo; 

■ Un nodo tiene siempre dos primos cercanos en cada dirección. Si no existe, debe ser 
creado como hoja virtual. 

■ Un nodo puede ser quitado sólo si son quitados todos sus hermanos y sólo si no es 
el primo cercano de un nodo existente. 




Figura 4.1: Estructura de datos tipo árbol graduado dinámico unidimensional. 



4.1.9. Análisis del error 



El error global entre los valores puntuales de la solución exacta en el nivel L, u^ x , y los 
valores de la solución numérica por multiresolución con un nivel máximo L, uj¿ R , puede 
ser descompuesto en dos errores [26, 32]: 

\\ u ex ~ U MR\\ ^ \\ u ex ~ u Fv\\ + \\ U FV ~ M Aíñll> (4.28) 

donde || • || es la norma £ , C 2 , o C°° . El primer error del lado derecho de (4.28), llamado 
error de discretización, es el error del esquema de volúmenes finitos en malla fina, para 
un nivel máximo L. Puede ser acotado por 

11*4 - u fv\\ < C2-^ L , O 0, (4.29) 

donde £ es el orden de convergencia del esquema de volúmenes finitos. En este caso, se 
utilizarán esquemas de segundo orden (en tiempo y espacio). Luego £ = 2. 

El segundo error del lado derecho de (4.28) es llamado error de perturbación. En [14] 
se prueba que si los detalles en un nivel de multiresolución k son truncados bajo cierta 
tolerancia prescrita £k, si el operador de evolución temporal discreto E es contractivo en 
la norma correspondiente, y si la tolerancia prescrita en el nivel k es 

e k = 2( fc " ¿ ) e , 
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entonces la diferencia entre la solución por volúmenes finitos en la malla fina y la solución 
obtenida mediante un algoritmo de multiresolución, se acumula en el tiempo y verifica 

\Wmr - UfvW < Cne, C > 0, (4.30) 
donde n es el número de pasos temporales. Considerando un tiempo fijo T = nAt, esto es 

14-^1^^, c>o. 

Para la ecuación lineal de convección-difusión, de (4.27), el paso temporal Ai debe verificar 

Ax 2 

At ^ 



4v + cAx 



Si se denota por X a la longitud del dominio, Ax al paso espacial en la malla fina, y en el 
caso de que la raíz del árbol graduado contenga sólo un nodo, se tiene Ax = X2~ L . Luego 

At = c J^ = C j , „ <C <1. 

4v + cAx 4v + cX2~ L 
Si se quiere que el error de perturbación sea del mismo orden que el error de discretización, 

e/At oc 2^ L . 

Por lo tanto, 

e2 2L (4v + cXT L ) oc 2~t L , 
y si se define el número de Peclet como Pe = ~, 

2"(í+ 1 ) ¿ 

e oc , r , . (4.31) 

Pe + 2( L + 2 ) v ' 

Para el caso invíscido {Pe — * oo), 4.31 es equivalente a los resultados obtenidos en [14]: 

e oc 2-^ L . 

Con esto, elegiremos una tolerancia de referencia: 

2-tt+i)£ 

e R = C -r—r. (4.32) 

Pe + 2( L + 2 ) V ' 

4.1.10. Cálculo del flujo conservativo 

Considérese una hoja Qk+i,2j+l con primos virtuales Qk+i,2j+2 y ^fc+i,2j+3 a la derecha. 
Su padre fi&j+i es una hoja. Como se ve en la figura 4.2, el flujo que sale de O^+i^j+i hacia 
la derecha i 7 A ; +i,2j+i^fc+i,2j+2 n o está en equilibrio con el flujo que sale de &k,j+i hacia 
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2/tí 


2j+3 






J 


1*1 





k 



Figura 4.2: Flujo saliente y entrante para dos niveles diferentes. 

la izquierda Fj-j+i-^kj- Es posible calcular directamente los flujos que salen de &k+i,2j+i 
hacia fijfc¿+i o pueden calcularse sólo los flujos en el nivel k + 1 y para determinar el flujo 
entrante a la hoja en el nivel k, éste será igual a la suma de los flujos salientes de las hojas 
en el nivel k + 1. 

Esta elección asegura una conservatividad estricta en el cálculo de los flujos entre 
volúmenes de control de niveles diferentes, sin un aumento significativo de las evaluaciones 
(generalmente costosas) de los flujos. 

4.1.11. Implementación del algoritmo 

A continuación se presenta la estrategia a seguir por el algoritmo. En primer lugar, 
dependiendo de la condición inicial dada, se crea un árbol graduado inicial. Luego se 
realiza la evolución temporal sobre las hojas y finalmente se actualiza el árbol graduado. 

■ Inicialización de parámetros: tiempo de simulación, tamaño del dominio, niveles de multire- 
solución, número de puntos en la malla fina, condición CFL, etc. 

■ Creación de la estructura de árbol graduado inicial: Cálculo de detalles mediante transfor- 
mada de multiresolución, obtención de la representación puntual csparsa. 

■ Evolución temporal: Cálculo del operador discreto de divergencia en todas las hojas, cálculo 
de la evolución temporal Runge-Kutta. 

■ Si algún valor resulta overflow, el proceso se considera numéricamente inestable. 

■ Actualización de la estructura de árbol. 

■ Estudio de distintos indicadores de error. Cálculo de la tasa de compresión. 

Notar que el algoritmo puede resumirse esquemáticamente por 

u n+1 = É(At) • rVT 1 • Tr(e) • M • u n , (4.33) 

donde M es el operador de multiresolución (Codificación), Tr(e) es el operador de trunca- 
miento con la tolerancia prescrita s, y E(Aí) es el operador discreto de evolución temporal. 
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4.2. Resultados numéricos 



En esta sección se reproducen los resultados numéricos en ID obtenido por Roussell 
et al. [32] y Bihari [ ] , utilizando para la evolución temporal un método explícito Runge- 
Kutta TVD de segundo orden; para la discretización del término advectivo se utiliza 
el esquema clásico de Roe, con los estados de izquierda y derecha obtenidos mediante 
interpolación ENO de segundo orden y para la discretización de la parte difusiva, se 
utiliza un esquema centrado de segundo orden. Se utiliza un orden de precisión para el 
interpolador de multiresolución de r = 2. Se utilizan mallas finas de 256, 512, 1024, 2048 
y 4096 volúmenes de control, tolerancias prescritas de e = 5 x 10 -3 y e = 10~ 3 , niveles de 
multiresolución hasta L = 13 y una estrategia para el operador de truncamiento Ek = ^rhk ; 
1 ^ k ^ L. 



4.2.1. Ecuación de convección-difusión en ID 

En el caso de que el flujo sea lineal, se considera la ecuación lineal de convección-difusión 
para (x, t) £ [—1, 1] x [0, oo[, c > 0, v > 0, 



du du d 2 
dt dx da 



(4.34) 



Si se considera como escala espacial característica al largo del dominio X y como escala 
temporal característica a T = c/X, (4.34) puede escribirse en la forma adimensional 
siguiente 

du du 1 d 2 u 
dt dx Pe dx 2 ' 

donde Pe denota el número de Peclet Pe = cX/u. Se estudia (4.35) asociada a la condición 
inicial 



u(x, 0) = uq(x) 



1, si x ^ 0, 
0, si x > 



(4.36) 



y condiciones de Dirichlet en la frontera 

u(-M) 

u(l,t) 

La solución analítica está dada por Hirsch [3 



1, 
0. 



u ex (x,t) = ^erfc (— 2~Y~T 



(4.37) 
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Se testearon tres casos en que el parámetro de control es el número de Peclet Pe: 

i) Pe = 100. En la figura 4.3 (izquierda) se muestra la solución numérica de (4.35) en el 
tiempo t = 0,3125. Se observa el fenómeno de propagación lineal de la discontinuidad 
hacia la derecha. Notar de la tabla 4.1, que los errores al comparar la solución 
obtenida mediante multiresolución y la solución obtenida sin aplicar el proceso de 
multiresolución, son bastante pequeños, pero se acumulan con el paso del tiempo. 

ii) Pe = 1000. En la figura 4.4 (izquierda) se muestra la solución numérica de (4.35) en 
el tiempo t = 0,5. La suavidad de la solución se debe principalmente a la difusividad. 

iii) Pe = 10000. Este caso es cercano al caso límite en que la viscosidad es baja en 
extremo, y el efecto "suavizante" es bastante lento. Este caso (y se verá lo mismo 
para el caso no lineal), es un ejemplo de que la solución invíscida puede obtenerse 
haciendo v — > 0. En la figura 4.5 (izquierda) se muestra la solución numérica de 
(4.35) en el tiempo t = 0,7031. Notar de la tabla 4.1, la tasa de compresión es 
considerablemente alta. 

Ecuación de convección-difusión, t=0.31 Coeficientes de ondelett significativos 




Figura 4.3: Izquierda: Solución inicial (rayas), solución analítica (linea), y solución numérica de 
multiresolución (círculos) en el tiempo t — 0,31 para la ec. de convección-difusión en ID asociada 
a la condición inicial (4.36), con Pe — 100, L — 7, e = 10~ 3 y Nq = 257. Derecha: Estructura de 
coeficientes de ondclcttc significativos correspondientes. 



42 



4.2 Resultados numéricos 



43 



Ecuación de convección-difusión, t=0.50 



1 0000000000 



0.8 



0.6 



0.4 



0.2 ■ 



-0.5 



-\®o- 



0.5 



o- O 



Coeficientes de ondelett significativos 



'§ 5 



a> 3 



-0.5 O 0.5 

x 



Figura 4.4: Izquierda: Solución inicial (rayas), solución analítica (linea), y solución numérica de multiresolución 
(círculos) en el tiempo í = 0,50 para la ec. de convección-difusión en ID asociada a la condición inicial (4.36), con 
Pe = 1000, L = 7, e = 10 -3 y Nq = 257. Derecha: Estructura de coeficientes de ondelette significativos. 



Ecuación de convección-difusión, t=0.70 
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0.5 



> 3 - 
z 



Coeficientes de ondelett significativos 
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Figura 4.5: Izquierda: Solución inicial (rayas), solución analítica (linea), y solución numérica de multiresolución 
(círculos) en el tiempo í = 0,70 para la ec. de convección-difusión en ID asociada a la condición inicial (4.36), con 
Pe = 10000, L = 7, e = 10 — 3 y No = 257. Derecha: Estructura de coeficientes de ondelette significativos. 
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Pe 


n 


A» 


ei 




&2 








1 nn 
1UU 


1U 


z4.uyuo 


íí nn vin- 
o.UU X 1U 


-i 


1 QQvIfl- 

l.ooX 1U 


-4 


n ni vi n- 
y.ui x iu 


-4 




100 


23.9254 


1.90x10" 


-i 


1.02x10" 


-4 


6.14x10" 


-4 




200 


23.6491 


4.31x10" 


-4 


6.58x10" 


-4 


8.65x10" 


-4 




600 


24.1358 


8.29x10" 


-4 


7.17x10" 


-4 


9.61 xlO" 4 


(*) 


1000 


10 


28.2134 


7.56x10" 


-1 


2.02x10" 


-4 


7.90x10" 


-3 




100 


27.6779 


8.38x10" 


-6 


6.80x10" 


-5 


6.56xl0~ 4 


(*) 




200 


28.7502 


1.45x10" 


-5 


2.72x10" 


-4 


9.77x10" 


-4 




600 


28.0683 


4.04x10" 


-1 


5.65x10" 


-4 


1.00x10" 


-3 


10000 


10 


32.0937 


1.77x10" 


-6 


2.37x10" 


-5 


5.82x10" 


-5 




100 


32.0901 


1.93x10" 


-5 


2.70x10" 


-4 


2.22x10" 


-4 




200 


32.0949 


1.82x10" 


-1 


5.72x10" 


-4 


4.43xl0~ 4 


(*) 




600 


32.1005 


2.94x10" 


-1 


7.16x10" 


-4 


9.23x10" 


-4 



Cuadro 4.1: Solución numérica de la Ecuación de Convección-difusión en ID, con condición 
inicial (4.36), L = 7, e = 10~ 3 y Nq = 257. Se adjuntaron figuras para los casos marcados 
con (*). 
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4.2.2. Ecuación de Burgers viscosa en ID 



Se llevaron a cabo experimentos con la ecuación de Burgers viscosa, la que contiene 
un término convectivo no lineal, para la cual se conoce solución analítica. Para (x, t) € 
[— 1, 1] x [0, oo[, la ecuación puede ser escrita en su forma adimensional: 

du d í v?\ 1 d 2 u 
dt dx \ 

donde Re = 4r es e l número de Reynolds. 



8t + dx\2) Redx 2 ' (4 " 38) 



Dato inicial suave 



Asociada a la ecuación (4.38), considérese la condición inicial 

u(x,0) = uo(x) = sin(7rx), — 1 ^ x < 1 (4.39) 

y condiciones de borde periódicas. Excepto en el caso límite cuando Re es muy grande, 
nunca existe un choque completamente discontinuo. Como se muestra en los resultados 
siguientes, se obtienen tasas de compresión cercanas a 4. 

Se presentan resultados para Re = 0,001, Re = 1, y Re = 10: 



i) Re = 0,001. Este caso corresponde a una difusividad grande, lo que provoca que el 
dato inicial se mantenga suave para todo tiempo t. Ver resultados en la tabla 4.2 y 
figura 4.6. 

ii) Re = 1. En n = 600 se advierte la creación de una N-onda y el diagrama de coefi- 
cientes de multiresolución es similar al obtenido en el caso invíscido. Ver resultados 
en la tabla 4.2 y figura 4.7. 

iii) Re = 10. Este caso produce resultados similares a los obtenidos en el caso invíscido. 
Debido a la capacidad del algoritmo de mantener perfiles afilados, la tasa de com- 
presión se mantiene bastante alta. La figura 4.8 muestra que el choque se encuentra 
en un estado casi estacionario. Ver resultados en la tabla 4.2 y figura 4.8. 



Notar que en todos los casos, los errores son bastante pequeños; por lo tanto la calidad 
de la solución no se ve comprometida al aplicar el proceso de multiresolución. 

En la sección siguiente se verá que para un número de Reynolds bastante grande, el 
problema viscoso no necesita un tratamiento especial, y puede utilizarse el proceso de 
multiresolución desarrollado para leyes de conservación hiperbólicas. 
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Re n 




ei 




e 2 








10 10 


3.0963 


1.54x10" 


-4 


2.24x10- 


-4 


5.12x10" 


-4 


100 


3.8254 


5.19x10" 


-4 


6.69x10- 


-4 


1.01x10" 


-3 


200 


4.6491 


6.23x10" 


-4 


6.72x10- 


-4 


7.49xl0" 4 


(*) 


600 


5.1358 


7.47x10" 


-4 


6.56x10- 


-4 


6.05x10" 


-4 


1000 


5.1358 


8.17x10" 


-4 


5.71x10" 


-4 


2.91x10" 


-3 


1 10 


4.0198 


9.86x10" 


-5 


1.62x10" 


-4 


4.21x10" 


-3 


100 


3.9876 


2.31x10" 


-4 


2.68x10" 


-4 


1.61x10" 


-4 


200 


3.9902 


2.92x10" 


-4 


2.74x10" 


-4 


3.48x10" 


-5 


600 


4.0299 


3.47x10" 


-4 


2.64x10" 


-4 


5.55xl0" 5 


(*) 


1000 


4.3742 


3.71x10" 


-4 


2.48x10" 


-4 


1.06x10" 


-4 


0.001 10 


4.0279 


1.15x10" 


-5 


1.84x10" 


-5 


4.52x10" 


-5 


100 


4.0198 


5.71x10" 


-5 


6.19x10" 


-5 


3.30x10" 


-4 


200 


4.0198 


7.42x10" 


-5 


1.67x10" 


-4 


4.79x10" 


-4 


600 


4.0021 


1.24x10" 


-1 


2.26x10" 


-4 


6.85x10" 


-4 


1000 


4.0021 


4.07x10" 


-4 


4.71x10" 


-4 


9.02xl0" 4 


(*) 



Cuadro 4.2: Solución numérica de la Ecuación de Burgers viscosa en ID, condición inicial 
(4.39), L = 7, e = 10" 3 y iV = 257. Se adjuntan fi guras para los casos marcados con (*). 
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Solución en n=1000, Re=0.001 , ecuación de Burgers viscosa Coeficientes de ondelett significativos 




x x 

Fi gura 4.6: Izquierda: Solución (rayas) y solución numérica de multiresolución (asteriscos) en el paso temporal 
n = 1000 para la ec. de Burgers viscosa, con Re = 0,001, L = 7, Nq = 257 y e = 10~ 3 . Derecha: Estructura de 
coeficientes de ondelette significativos correspondientes. 



Solución en n=600, Re=1 , ecuación de Burgers viscosa Coeficientes de ondelett significativos 




Fi gura 4.7: Izquierda: Solución (rayas) y solución numérica de multiresolución (asteriscos) en el paso temporal 
n = 600 para la ec. de Burgers viscosa, con Re = 1, L = 7, ./Vo = 257 y e = 10 -3 . Derecha: Estructura de coeficientes 
de ondelette significativos correspondientes. 
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Solución en n=200, Re=10, ecuación de Burgers viscosa 



Coeficientes de ondeleít significativos 



Figura 4.8: Izquierda: Solución (rayas) y solución numérica de multiresolución (asteriscos) en el 
paso temporal n — 200 para la ec. de Burgers viscosa, con Re = 10, L = 7, No — 257 y e — 10 -3 . 
Derecha: Estructura de coeficientes de ondelettc significativos correspondientes. 



Dato inicial discontinuo 



Asociada a la ecuación (4.38), considérese la condición inicial 

u(x, 0) = Uq(x) = 



1, si x ^ 0, 
0, si x > 



y condiciones de Dirichlet en la frontera 

u(-M) 

u{l,t) 

La solución analítica está dada por (4.22) 



1, 
0. 



1 



(4.40) 



(4.41) 



La solución numérica de (4.38) en el tiempo t = 0,5 se muestra en la parte izquierda de la 
figura 4.9 para Re = 1000, e = 10~ 3 y L = 7 escalas de multiresolución, correspondientes 
a un máximo de 512 volúmenes de control en la malla fina. En la parte derecha de la 
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figura 4.9 se representan los coeficientes de ondelette significativos. Es posible notar el 
efecto de una propagación no lineal del choque hacia la derecha, además puede notarse la 
difusividad cerca de la discontinuidad. 



Evolución Runge Kutta de orden 2, CFL=0.5, T=0.50 
* * * t-N»* 



Coeficientes de ondelette significativos 



Sol por MR 
- Sol exacta 
Sol inicial 




Figura 4.9: Izquierda: Condición inicial (rayas), solución analítica (linea), y solución con multi- 
resolución (asteriscos) en el tiempo t = 0,5, ec. de Burgers viscosa, Re — 1000, L — 7, N = 257 y 
e = 10~ 3 . Derecha: Estructura de coeficientes de ondelette significativos. 



También se presenta la solución numérica obtenida mediante el esquema ENO de se- 
gundo orden, con Runge-Kutta de segundo orden (EN02-RK2) pero sin aplicar multire- 
solución (parte izquierda de la figura 4.10). La evolución temporal de los errores entre las 
soluciones analítica y calculada mediante volúmenes finitos con y sin multiresolución se 
presenta en la parte derecha de la figura 4.10. Notar que los errores están bajo la tolerancia 
prescrita de e = 10~ 3 . Como una medida de la mejora en velocidad alcanzada mediante 
la utilización del análisis de multiresolución, se utiliza la tasa de compresión definida por 

" -Ñ¡¡$Wv (4 - 42) 

donde D n es el conjunto de coeficientes de ondelette significativos, en todos los niveles de 
multiresolución, en el paso temporal n. 

En las tablas 4.2.2 y 4.2.2 se muestra para diferentes tiempos la constante de proporción 
V entre el tiempo de CPU total para calcular la solución numérica sin multiresolución y el 
tiempo de CPU total para calcular la solución numérica con multiresolución. Nótese que 
de los resultados de las tablas se concluye que la solución numérica tarda alrededor de 1.6 
veces el tiempo de CPU que la solución de multiresolución. 
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Solución numérica sin mulliresolución 



O VF 
exacta 



0.5 



1 r 

0.9 - 

0.8 - 

0.7- 

0.6 - 

oT 0.5 ■ 

0.4 - 

0.3- 

0.2 - 

0.1 - 

- 




x ^ o" 3 Emires en norma || ■ ||i, t=0.5 



- VF 

- Multiresolucion 



0.1 0.2 0.3 0.4 



0.5 



Fi gura 4.10: Izquierda: Solución analítica (linea), y solución numérica sin multiresolucion (circuios) en el tiempo 
í = 0,5 para la ec. de Burgers viscosa, con Re = 1000, L = 7, Nq = 257 y e = 10 — 3 . Derecha: Errores entre las 
soluciones analítica y de volúmenes finitos con y sin multiresolucion. 
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0.3 


10.1484 




0.3 


1 Q858 




0.3 


9.8646 




0.4 


9.9543 




0.4 


10 6509 




0.4 


9.8646 




0.5 


9 39fi0 




0.5 


10 001 3 




0.5 


9.4993 


« 
o 


001 Q 


32 0^86 


q 


001 9 


99 31 96 


i n 


001 9 


98 1 088 

ZjO . 1UOO 




n i 

yj . 






n i 


1 691 7 




n 1 






0.2 


9.4997 




0.2 


9.2067 




0.2 


8.6282 




0.3 


9.3264 




0.3 


8.3413 




0.3 


7.0034 




0.4 


9.1604 




0.4 


7.5994 




0.4 


6.9844 




0.5 


9.0326 




0.5 


6.8858 




0.5 


6.3925 


11 


0.0019 


28.3031 


12 


0.0019 


28.3012 


13 


0.0019 


28.1505 




0.1 


8.2575 




0.1 


8.1348 




0.1 


8.1003 




0.2 


7.5302 




0.2 


7.4280 




0.2 


7.2222 




0.3 


7.3771 




0.3 


7.3220 




0.3 


7.1219 




0.4 


6.9207 




0.4 


6.6343 




0.4 


6.7472 




0.5 


6.2466 




0.5 


6.1026 




0.5 


6.1049 



Cuadro 4.3: Tasa de compresión para distintos niveles de multiresolución, hasta t = 0,5 
para la ecuación de Burgers viscosa en ID, condición inicial (4.40). 
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t 


V 


0.06 


2.0011 


0.12 


1.9912 


0.18 


1.8123 


0.24 


1.7780 


0.36 


1.6761 


0.42 


1.6302 


0.48 


1.6079 



Cuadro 4.4: Proporción V entre el tiempo de CPU total de la solución numérica EN02 en malla 
fina y el tiempo de la solución de multiresolución. N = 257, L = 7 y e = 10~ 3 . 

Al aumentar el número de puntos en la malla fina, los resultados obtenidos son aún 
mejores, y en este caso la solución de multiresolución tarda menos de la mitad del tiempo 
total de CPU que tarda la solución numérica que no utiliza multiresolución. 



t 


V 


0.06 


3.0444 


0.12 


2.6358 


0.18 


2.5129 


0.24 


2.5089 


0.36 


2.4761 


0.42 


2.4341 


0.48 


2.4192 



Cuadro 4.5: Proporción V entre el tiempo de CPU total de la solución numérica EN02 en malla 
fina y el tiempo de la solución de multiresolución. N = 513, L = 9 y e = 10~ 3 . 

En el caso de sistemas de leyes de conservación o en el caso de problemas multidimen- 
sionales, se espera que V sea aún más significativo. 
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4.2.3. Ecuación de reacción-difusión en ID 



Otro prototipo de una ecuación parabólica no lineal es la ecuación de reacción- difusión. 
En este caso, la no linealidad no está más en el término advectivo (como en la ecuación 
de Burgers viscosa) sino en el término fuente. Para (x,t) £ [0,20] x [0,oo[, la ecuación 
puede ser escrita en su forma adimensional: 

du d 2 u . . A 1 . 

* = b? + w< (4 - 43) 

con 

S(u) = ^(l-n)exp ^ (1 ~ U) 1 , (4.44) 
2 a{± — u) — 1 

donde a es la tasa de temperatura y (3 es la energía de activación adimensional (número 
de Zeldovich). Se estudia (4.43) asociada a la condición inicial 



u(x, 0) = uq(x) = < 



1, si x ^ 1, 

(4.45) 

exp(l — x), si x > 1. 



Esta ecuación conduce al modelo de la propagación de una llama premezclada en ID, 
donde las difusividades de masa y calor son iguales. La función u representa la temperatura 
adimensional, que varía entre y 1. La masa parcial de gas sin quemar es 1 — u. Se elige 
una condición de Neumann en la frontera izquierda y una condición de Dirichlet en la 
frontera derecha. 

du . 

m M = o, 

u(20,í) = 0. 

Los parámetros son a = 0,8 y (3 = 10. EL tiempo final (adimensional) es í/ = 10. En este 
ejemplo, la no linealidad del término fuente implica que Ai ~ O(Ax). 

La velocidad de la llama, definida por 

v f = í Sdx (4.46) 
Jn 

se compara con los valores asintóticos dados por Peters 8¿ Warnatz [32]. 

En la figura 4.12 se observa la propagación de la llama en la dirección x. El mayor nivel 
es alcanzado en la región de la zona de reacción, es decir, para x w 10. 
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Propagación de llama premezclada, CFL=0.1, 
1 .8 




Coeficientes de ondelette significativos, No=513 

7 i 1 m 1 1 



Figura 4.12: Izquierda: Condición inicial (rayas) y S(u) inicial (puntos), solución numérica sin 
multiresolución (linea), solución numérica con multiresolución (asteriscos) y S(u) (puntos-rayas), 
en el tiempo t = 10 para la ec. de reacción-difusión, con a = 0,8, (3 = 10, L = 7, Nq = 513 y 



- 10~ 3 



Derecha: Estructura de coeficientes de ondelette significativos, t = 0,5. 



Método 


v f 




VF 


0.9146 




MR e = 5 x 1(T 2 


0.9182 


12.5648 


MR e = 1(T 3 


0.9151 


13.8977 (*) 


Valor asintótico 


0.9080 





Cuadro 4.6: Velocidad de la llama y tasa de compresión para la solución numérica de 
(4.43) sin multiresolución (VF), y a dos niveles distintos de tolerancia prescrita para el 
caso multiresolutivo. Nq = 513. (*) representado en la figura 4.12. 
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Capítulo 5 



Ecuación de convección-difusión 
fuertemente degenerada 

En este capítulo se presentará un método numérico para obtener soluciones aproxi- 
madas de problemas provenientes de la sedimentación de suspensiones floculadas. Estos 
procesos se utilizan para lograr la separación de una suspensión de pequeñas partículas 
suspendidas en un líquido viscoso, en sus componentes sólido y líquido bajo la acción de 
la fuerza de gravedad. Estos procesos se usan ampliamente en la industria minera, por 
ejemplo para recuperar el agua de las suspensiones que salen de los procesos de flotación 
[11]. 

La idea principal es aplicar los métodos de multiresolución a los esquemas desarrollados 
por Bürger et al. [5, 7, 8, 9, 10] y observar que el método de multiresolución descrito y ejem- 
plificado en los capítulos anteriores es de gran ayuda para reducir el costo computacional 
en este tipo de problemas sin afectar la calidad de la solución. 

Se dará una breve descripción del problema físico y su modelación mediante una ley de 
conservación fuertemente degenerada con flujo no lineal [7]. El efecto de la compresibilidad 
del sedimento puede ser descrito por un término difusivo fuertemente degenerado, mien- 
tras el flujo unidimensional contribuye una discontinuidad de flujo a la ecuación parcial 
diferencial. Se presentará un esquema de segundo orden desarrollado en Bürger y Karlsen 
] para resolver este tipo de problemas y finalmente se desarrollan ejemplos numéricos 
para comparar con los resultados publicados en [7, 8, 9]. 

Considérese el caso de una suspensión floculada en un ICT (Ideal Continuous Thi- 
ckener) como el de la figura 5.1, derecha. Un ICT es un espesador cilindrico sin efectos 
de pared, en que las variables dependen sólo de la altura x y el tiempo t. En x = H se 
tiene una superficie de alimentación y en x = se tiene una superficie de descarga, lo que 
produce una operación continua del proceso. Esta modelación es prácticamente obsoleta, 
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pero es de gran utilidad al momento de ejemplificar el comportamiento simplificado de los 
procesos de sedimentación. El caso especial de sedimentación batch se muestra en la parte 
izquierda de la figura 5.1. El recipiente es cerrado. 



overflow 




Figura 5.1: Izquierda: Columna de sedimentación Batch. Derecha: ICT (Ideal Continuous Thi- 
ckener) [8]. 

En el caso unidimensional, la teoría de la sedimentación produce ecuaciones de equili- 
brio de masa y momentum lineal que pueden simplificarse [11] hasta obtener una ecuación 
parabólica fuertemente degenerada de la forma 

d t u + d x f{u) = dl x A{u), (5.1) 

con (x,t) e]0,l[x[0,T[ y el coeficiente de difusión integrado dado por 

A(u) := / a(s)ds, a(u) ^ 0. (5.2) 
Jo 

En general, se permite que el coeficiente de difusión a{u) sea cero sobre intervalos de u. En 
tales casos, (5.1) es una ecuación hiperbólica. Por esto, (5.1) se denomina también ecuación 
hiperbólica-parabólica. Aún cuando este tipo de ecuaciones modelan una gran variedad de 
fenómenos, se enfatizará en las aplicaciones a los procesos de sedimentación-consolidación. 

Las soluciones de (5.1) desarrollan discontinuidades debido a la no linealidad de la 
función de densidad de flujo f(u) y a la degeneración del coeficiente de difusión. Esto 
lleva a considerar soluciones entrópicas para tener un problema bien puesto. Aún más, 
cuando (5.1) es puramente hiperbólica, los valores de la solución se propagan sobre rectas 
características que podrían intersectar las fronteras del dominio espacio-tiempo desde el 
interior, y esto requiere tratar a las condiciones de Dirichlet como condiciones entrópicas 

[']■ 
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Una gran parte de las ecuaciones constitutivas que se proponen para estos procesos, 
implican que a{u) tiene un comportamiento degenerado, es decir, a(u) = para u ^ u c 
y a{u) salta en u c a un valor positivo, donde u c es una constante llamada concentración 
crítica. Se enfatiza entonces el hecho de que el coeficiente de difusión a{u) es degenerado, 
lo que hace evidente la naturaleza hiperbólica-parabólica de la ecuación diferencial (5.1). 



Considérese el problema de valores iniciales y de frontera (PVIF) siguiente 



d t u + d x {q{t)u + /(«)) 
u(x, 0) 
u(H, t) 

f(u(0,t)) - d x A(u(0,t)) 



d 2 xx A(u), (x,t) €]0,H[x[0,T[, 
u (x), x £ [0,H], 
0, ÍG]0,T] 

o, íe]0,r], 



conocido como el Problema A. Considérese además el Problema B 

a2 



d t u + d x (q(t)u + f(u)) 
u(x, 0) 

q{t)u{H,t)-d x A{u{H,t)) 
f(u(0,t)) - d x A(u(0,t)) 



d z xx A{u), (x,t)e)0,H[x[0,T[, 
u (x), x £ [0,H], 
*(í), te]0,T] 

o, íe]0,r]. 



(5.3) 
(5.4) 
(5.5) 
(5.6) 



(5.7) 
(5.8) 
(5.9) 
(5.10) 



Para ambos problemas, / se supone continua y diferenciable a trozos, / ^ 0, sop (/) C 
[0,u m áx], H/'Hoo < oo, a(u) ^ 0, sop (a) C sop (/), a(u) = para u ^ u c , < u c < u máx , 
q{t) ^ 0, Vi G [0,T], !TV(g) < oo, IV(g') < oo. 

En [7] se prueba la existencia y unicidad de solución entrópica para cada uno de estos 
problemas. 

En los modelos de sedimentación-consolidación de suspensiones floculadas, la coorde- 
nada x aumenta verticalmente, u = u{x, t) representa la concentración volumétrica sólida 
local, q(t) ^ es la velocidad media del flujo de la mezcla (puede ser controlada exter- 
namente), f(u) es una función dada que relaciona la velocidad relativa local sólido-fluido 
con la concentración de sólidos local, y 



a[u) 



f{u)a' e {u) 
Aggu 



(5.11) 



donde Ag > denota la diferencia de densidad de masa sólido-fluido, g es la aceleración 
de gravedad, y o~' e (u) ^ es la derivada de la función de rigidez sólida efectiva. 

La propiedad de mayor interés, es que generalmente se supone el siguiente comporta- 
miento para a e {u): 



a e (u) < 



= cte., si u ^ u c , da e 
> U, SI u > u c , 



= 0, si u ^ u c , 
> 0, si u > u c - 



(5.12) 
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Notar que la naturaleza degenerada de la ecuación diferencial (5.1) es heredada de esta 
propiedad. 

Las propiedades materiales específicas de la suspensión son descritas por f(u) y a e (u). 
Ejemplos típicos para estas funciones modelo son la función de densidad de flujo del tipo 
Michaels and Bolger [8] 

c 



f(u) = VooU 1 



íí 



, Voo < 0, C > 1 



(5.13) 



y la función de rigidez sólida efectiva ley de potencia 

0, si u < u r 



a e (u) 



ero 



ao > 0, n > 1. 



(5.14) 



1 ) , si u > u c 



Las condiciones (5.4) corresponden a una distribución inicial de concentración dada, 
la condición (5.5) corresponde a prescribir el valor de la concentración en x = L, las 
condiciones (5.6) y (5.10) equivalen a reducir la densidad de flujo del volumen sólido en 
el fondo del recipiente a su parte convectiva q(t)u(0,t) y la condición (5.9) corresponde a 
una condición de flujo en x = L. 



5.1. Esquemas de segundo orden 

Para el esquema explícito a desarrollar, se utilizará una discretización similar a la 
utilizada en la sección 4.1.4 (ver detalles en [7]). Los términos advectivo y difusivo son 
aproximados de diferente forma, con el fin de obtener una discretización que mantenga la 
conservatividad en ambos términos. Para la parte advectiva puede utilizarse el esquema 
de Roe clásico con una interpolación ENO de segundo orden, ya utilizado en los capítulos 
anteriores, o bien puede utilizarse un esquema de Engquist-Osher [1 8] modificado para ser 
de segundo orden [8, 9, 19]. Para la parte difusiva, se necesita un esquema centrado de 
segundo orden que mantenga la conservatividad [ ]. 

Dado que el principal interés se encuentra en la discretización del término difusivo, 
considérese la siguiente ecuación puramente difusiva: 

dtu = d1 x A{u), (5.15) 



A(u) = / a(s)ds. (5.16) 
Jo 

Una formulación conservativa de diferencias finitas para esta ecuación es 



n+l 

Ai (Ax) 



ir; • ' - ,<■; A(^_J - 2A(u]) + A(u] +1 ) 



(5.17) 
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Este esquema es estable y convergente bajo la condición CFL (ver [9]) 



2max|a(u)| /A * ^ 1. (5.18) 
u (Ax) ¿ 



Además, debe recordarse que el esquema explícito utilizado para la ecuación puramente 
hiperbólica es estable bajo la condición CFL (ver [21]) 

máx|/'(n)|^ ^ 1. (5.19) 
u Ax 

El esquema interior resultante para la ecuación (5.1) (si se utiliza el esquema de Roe 
clásico (4.12) con una interpolación ENO de segundo orden para la parte advectiva) es: 



— — t + q(nAt) J 2 — = . - ^ J - J -- 

At Hy ' Ax Ax (Ax) 2 



con 



F 



2 



La evolución temporal se hará mediante el método Runge-Kutta de segundo orden 
utilizado en los capítulos anteriores. 

Las condiciones de borde (5.6) y (5.10) prescritas en x = se discretizan utilizando 
(5.20) haciendo: 

f(u(0, 0) - d x A(u(0, t n )) m F n i — = Q) (5 22) 

~2 Ax 

de donde se obtiene la expresión para la actualización del flujo en Uq 

< +1 - < + „(„A t )^ + ? = A W» - . (5.23) 



Ai v ' Ax Ax (Ax) 

Esta formulación evita utilizar un valor artificial u n L 1 . 

Para el problema A, la condición de borde en x = H se aproxima simplemente poniendo 
u^r Q = 0, en cambio para el problema B, (5.9) se aproxima haciendo 

Akñr ,,)-A(u n N ) 
q(nAt)u% + F- Q+h - No+1> Ax No> = *(nAí). (5.24) 

Con esto, se obtiene la expresión para la actualización del flujo en 

n+l_.,n ,t,íV, A+\ _ „( m K+\n,n F 



U N U N 



, *(nAt)-q(nAt)u n No ^ -í _ Ajvfr^) - A(u\) 



At Ax Ax (Ax 
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Como alternativa a la discretización de la parte advectiva, puede utilizarse un esquema 
de Engquist-Osher modificado mediante extrapolación de variables MUSCL (Monotonic 
Upwind Scheme for Conservation Laws) para lograr un esquema de segundo orden [8, 9, 
19, 21, 20]. Para ello se introduce una función u n {x) lineal a trozos definida por 

U n (x) = U] + s](x - Xj), X e]Xj_i/2, x j+1/2 [, 

donde s™ es una pendiente adecuada, construida a partir de u n . En las regiones donde 
Sj = 1, la reconstrucción es lineal y el error de truncamiento es 0((Ax) 2 ). En las regiones 
donde s™ = 0, la reconstrucción es constante a trozos y el error de truncamiento es 
O(Ax). Es necesario utilizar limitadores de pendiente para forzar la monotonía de la 
reconstrucción. En este caso, se utilizará el 8— limitador (ver [21, 22]) 

(n i Tí n i TI r¡ i TI i¡ i TI n , TX n i TI \ 

»t ■ 

donde MM es otra función tipo Min-Mod definida por 

mín(a, b, c), si a, b, c > 0, 
MM(a,b,c) := { máx(a,6,c), si a, b, c < 0, (5-26) 
0, e.o.c. 

Luego se extrapola la información hacia la frontera de cada volumen de control, con lo 
que 

uf:=«7-^-y, uf : <<; • A ^'¡- (5-27) 
Así, el esquema upwind interior de segundo orden correspondiente se escribe 

Ai + ^ nAt > + Ax " {Axj* ' (5 ' 28) 

donde f EO (Uj, Uj +1 ) := / + (u") + / _ (u™ +1 ) es el flujo numérico de Engquist-Osher [18], 
/+(«) = /(0)+ / máx(/'( S ),0)d S , /-(«)=/ mm{f , (s),0)ds. (5.29) 



o 



máx| /'(«)!— + 2máx|a(u)| 7x:: ^ < 1. (5.30) 



Este esquema es estable bajo la condición CFL (ver [19]) 

Ai , .... Ai 

h 2 max \a(u)\-— A — 

Ax u 1 y '\Ax 

Las condiciones de borde (5.6) y (5.10) prescritas en x = quedan entonces 

Ai + q[n¿Xt) ^T + Ax " (AxY ( j 

y la condición de borde (5.9) queda 

< l ~ "ftp *(n*t) - q(nAt)u" No f E °(u% _ lt u% ) _ Ajn^) - A(u« No ) 

At Ax Ax (Ax) 2 ' { ' 
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5.2. Un algoritmo de multiresolución 

Se presenta a continuación una breve descripción de un algoritmo de multiresolución 
para resolver numéricamente una ecuación parabólica fuertemente degenerada. 

1. Inicialización de parámetros y variables: 

■ Longitud del dominio H, 

■ concentración crítica u c , 

■ orden de la interpolación de multiresolución r, 

■ niveles de multiresolución L, 

■ número de puntos y paso en la malla fina Nq y ho, y en cada nivel N¡¿ y hk, 

■ tolerancia prescrita e y estrategia de truncamiento 

■ tiempo de simulación tf, 

■ constantes de Lipschitz para a(u) y f'(u), 

■ condición CFL: 

, .... ..Ai , . . ,.Aí 

max|j (u)\- h 2 max \a(u)\ ^ 1. 

■ paso temporal Ai, 

CFL • h 



At 



máxu |/'(íí)| + 2máx u \a(u)\/ho 

■ condiciones iniciales uq y 

■ otros parámetros del modelo (5.38): Voo, C, n, tt m á x , Ag, etc. 

■ Inicialización de la estructura de datos. (En este caso, estructura esparsa). 

2. Aplicación de la codificación a la condición inicial: Este proceso entrega los coefi- 
cientes de ondelette significativos y los valores de la solución en las posiciones corres- 
pondientes a coeficientes de ondelette significativos. Se incluyen los safety points. 

3. Evolución temporal: Se utiliza un método Runge-Kutta de segundo orden. 

■ Primer paso intermedio Runge-Kutta, 

■ Segundo paso Runge-Kutta, 

■ Actualización de los flujos y actualización de la solución, 

■ Imposición de condiciones de contorno fijas y condiciones de flujo, 

■ Se aplica el paso 2. a la solución actual y se itera hasta alcanzar el tiempo final. 

4. Salidas: Se realizan gráficos de la solución numérica y coeficientes de ondelette sig- 
nificativos correspondientes. Se calculan además tasas de compresión y tiempos de 
CPU para comparar con la resolución obtenida sin utilizar multiresolución. 
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5.3. Ejemplos numéricos 



Se calculan soluciones de los problemas A y B utilizando los esquemas numéricos des- 
critos en la sección anterior, con una discretización para el flujo de tipo Enqguist-Osher, 
dada por (5.28). Se reproducen algunos resultados numéricos obtenidos por Bürger et al. 
[7, 8, 9] y Bustos et al. [11]. 



5.3.1. Sedimentación batch de suspensión ideal 

Considerar en primer lugar, el proceso de sedimentación batch de suspensión ideal en 
una columna de asentamiento [15]. El caso ideal permite formular el proceso como 

du 



Ot 



+ 



df(u) 

dx 
u(x, 0) 

u(0,t) 
u(L,t) 



0, x G M, t > 0, 
uq(x), x £ [0,H[, 



UL, 



t > 0, 
t > 0. 



En el ejemplo se considera una columna de asentamiento de longitud H = 1, una con- 
centración inicial uo(x) = 0,25, condiciones de borde Uqo = 0,642 y no = 0. Se elige una 
ecuación constitutiva para la función de densidad de flujo sólido. Se utiliza la función 
descrita por Shannon (1963, consultar [11]) 

/(u) = (-0,33843u + l,37672tí 2 - l,62275u 3 - 0,11264u 4 + 0,902253w 5 ) x 1CT 2 [m/s]. (5.33) 




0.2 0.3 



Figura 5.2: Función de densidad de flujo f(u) para el problema de sedimentación batch de sus- 
pensión ideal. Unidad: [m/s]. 



En las figuras 5.3-5.5 se muestran soluciones numéricas para t = 60 [s], t = 300 [s] y t = 
3600 [s] obtenidas mediante el esquema de segundo orden descrito en la sección anterior, 
aplicando multiresolución. En t = 3600 [s] la solución ya alcanzó un estado estacionario. 
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Suspensión ideal, columna de asentamiento, É=60[s] Coeficientes de ondelette, IVo =257, L = 5. e = 10^ 




Figura 5.3: Izquierda: Condición inicial (rayas) y perfil de concentración a í = 60[s] para el problema de sedi- 
mentación batch de suspensión ideal (Asteriscos). Derecha: Coeficientes de ondelette significativos correspondientes. 



Suspensión ideal, columna de asentamiento, t =30Q[s] Coeficientes de ondelette, N =257, ¿ = 5, £ =10^ 




Figura 5.4: Izquierda: Condición inicial (rayas) y perfil de concentración a í = 300[s] para el problema de sedi- 
mentación batch de suspensión ideal (Asteriscos). Derecha: Coeficientes de ondelette significativos correspondientes. 
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Figura 5.5: Izquierda: Condición inicial (rayas) y perfil de concentración a t = 3600 [s] para 
el problema de sedimentación batch de suspensión ideal (Asteriscos). Derecha: Coeficientes de 
ondelette significativos correspondientes. 



En la tabla 5.1 se muestran la proporción V, tasa de compresión y errores entre la 
solución calculada utilizando multiresolución y la solución calculada sin multiresolución 
(ver sección 3.5). 



t[s] 


V 




ei 




&2 








60 


4.3457 


7.8456 


2.64x10" 


-5 


6.54x10" 


-6 


9.03x10" 


-6 


300 


5.6212 


5.8456 


1.70x10" 


-5 


6.39x10" 


-6 


1.12x10" 


-5 


1800 


5.9443 


14.9168 


7.28x10" 


-5 


2.98x10" 


-5 


4.35x10" 


-5 


3600 


6.1385 


29.8479 


8.89x10" 


-5 


4.04x10" 


-5 


6.50x10" 


-5 



Cuadro 5.1: Sedimentación de suspensión ideal, e = 1,0 X 10 4 , Nq — 257 y L = 5. 
Notar que los errores permanecen siempre bajo la tolerancia prescrita e = 1,0 x 10~ 4 . 
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5.3.2. Caso batch de suspensiones floculadas: primer ejemplo 

En este ejemplo se considera el caso batch de suspensión homogénea de concentración 
inicial uq(x) = 0,15 en un a columna de asentamiento cerrada, es decir, se considera el caso 
de q = 0, con una concentración prescrita en x = 1 dada por (5.5). El dominio espacial 
es [0, 1] y la concentración crítica es u c = 0,23. Notar que la discontinuidad entre u = 
y u = uq es un choque. Aún más, el problema (5.3)-(5.6) es un problema de Riemann, 
en el sentido de que el dato inicial consiste en dos estados constantes y la solución, en 
general, consistirá en ondas elementales: choques, ondas de rarefacción y discontinuidades 
de contacto [ ] . 

Como función de densidad de flujo, se utiliza una función Kynch batch Richardson-Zaki 
con parámetros correspondientes a suspensión de cobre [8]. 

f(u) = -6,05 x 10~ 4 u(l - u) 12 ' 59 [m/s]. (5.34) 

Se utilizará la función a' e (u) dada por ([9, 11]) 

a' e (u) = — (I00(n/n c ) 8 - l) [Pa], si u > u c . (5.35) 

Luego 

0, si u ^ u c = 0,23, 

f (¿)>a], siu>u c . 



a' e {u) = < 



(5.36) 



La función a{u) (5.11) está dada entonces por 

( 

a(u) = 

4,84X11 

u%Agg 

con Ag = 1500 [Kg/m 3 ] y g = 9,81 [Kgm/s 2 ]. 



0, si u ^ u c = 0,23, 

4,84xl0" 1 M 7 (l-íi) 12i59 



(5.37) 



La figura 5.6 muestra las funciones modelo f{u) y a{u). La función A(u) correspondiente 
al término difusivo integrado, se calcula mediante las fórmulas (5.39)-(5.40). 

En la tabla 5.2 se muestran la proporción V, tasa de compresión y errores entre la 
solución obtenida utilizando multiresolución y la solución obtenida sin multiresolución. 

Notar de la tabla 5.2, que los errores se encuentran por debajo de la tolerancia prescrita. 
Notar además los excelentes resultados en cuanto a proporción V (correspondiente al 
tiempo total de CPU en ambos casos). Los resultados en cuanto a tasa de compresión no 
son excelentes, pero hay que tomar en cuenta que se está considerando una malla de 129 
puntos. 
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O 0.2 0.4 0.6 0.8 1 O 0.2 0.4 0.6 0.8 1 

u u 

Figura 5.6: Funciones modelo f(u) (izquierda) y a(u) (derecha) para el problema de sedimentación- 
consolidación. Las unidades son [m/s] para f(u) y [m 2 /s] para a(u). 



t[s] 


V 




ei 




&2 




&00 




60 


6.5737 


17.8796 


1.29x10" 


-4 


8.72x10- 


-5 


5.33x10- 


-5 


1800 (*) 


5.7349 


9.4132 


1.99x10- 


-4 


9.06x10- 


-5 


7.42x10- 


-5 


3600 (*) 


6.1982 


9.1246 


2.77x10- 


-4 


2.67x10- 


-4 


9.61x10- 


-5 


7200 (*) 


6.2110 


9.1246 


3.21x10- 


-4 


4.67x10- 


-4 


2.41x10- 


-4 


14400(*) 


7.9244 


9.4132 


8.92x10- 


-4 


7.81x10- 


-4 


6.18x10- 


-4 



Cuadro 5.2: Suspensiones floculadas, primer ejemplo. Multiresolución utilizando e = 10 , 
iVo = 129 y L = 5. (*): figuras 5.7 - 5.9. 

En la figura 5.7 se presenta un perfil de concentración en un tiempo t = 1800[s], uti- 
lizando multiresolución. La solución se calcula utilizando 129 puntos en la malla fina, 
con una estrategia de truncamiento Ek = ^r-k- Se presenta además la configuración de 
los coeficientes de ondelette significativos. Notar que cuanto más perfilada es la discon- 
tinuidad, menor es el número de coeficientes de ondelette significativos asociados a tal 
discontinuidad. 

En la figura 5.8 se presenta un perfil de concentración en un tiempo í = l[h], utilizando 
multiresolución. La solución se calcula utilizando 129 puntos en la malla fina. Se presenta 
además la configuración de los coeficientes de ondelette significativos correspondientes. 
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Casa Bateh. sedimen 



ompresíon t =18QO[.s] 



É=0.5[/i] (MR) 
-t=0.5[h] (DF) 




t .'< .ii' Hc-ici utos de (.(íidokittí. 1 s i^nitíc ¡\ li\i>s , Aj, =129, £ = 



Ni-w-íles de Multiiosolucioi 



Figura 5.7: Izquierda: Condición inicial (rayas) y perfil de concentración ai = 1800[s] para el pro- 
blema de sedimentación-consolidación (asteriscos). Derecha: Coeficientes de ondclette significativos 
correspondientes. e=10~ 3 , AT = 129 y L = 5. 



I Bil.tfl). SOdílllO 



0. 8 
0.7 1 
0. 6 : 

L 0.5 




de ondelotto s ijfiiilic a 



, Ai, =129, ¿ = 5, e^lO" 3 



Ni w los do Muimos olut 



Figura 5.8: Izquierda: Condición inicial (rayas) y perfil de concentración a t = 3600 [s] para el pro- 
blema de sedimentación-consolidación (asteriscos). Derecha: Coeficientes de ondelette significativos 
correspondientes. e=10~ 3 , A r = 129yL = 5. 



En la figura 5.9 se presenta un perfil de concentración en un tiempo t = 4[/i], utilizando 
multiresolución. La solución se calcula utilizando 129 puntos en la malla fina. Se presenta 
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además la configuración de los coeficientes de ondelette significativos correspondientes. En 
este tiempo, la solución ya se encuentra en un estado estacionario (ver además 5.10 y [7]). 




£=14400[s] Coeficientes de ondelette s ignific atKos , Afc =129, L = 5, e=10^ 



1 

0.9 

0.8 

0.7 - 

0.6 - 

c 0.5 - 

0.4 - 

0.3 

0.2 

0.1 - 

— 
1 



3 4 5 

Ni\eles de Mtiltiies olucu 



Figura 5.9: Izquierda: Condición inicial (rayas) y perfil de concentración a, t — A[h] para el pro- 
blema de sedimentación-consolidación, caso Batch. (asteriscos). Derecha: Coeficientes de ondelette 
significativos correspondientes, e = 10~ 3 , Nq = 129 y L = 5. 




Figura 5.10: Perfiles de concentración hasta t — 12[h] para el problema de sedimentación- 
consolidación, caso Batch. e — 1CP 3 , N = 129 y L = 5. 



Los resultados numéricos concuerdan con los resultados obtenidos por Bürger et al. [7] . 
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5.3.3. Caso batch de suspensiones floculadas: segundo ejemplo 

En este ejemplo se considera el caso batch de suspensión homogénea de concentración 
inicial uq(x) = 0,05 en una columna de asentamiento cerrada {q = 0) de menor longitud: 
H = 0,16[m] (ver [9]). La concentración crítica es u c = 0,07. 

Como función de densidad de flujo, se utiliza la función f(u) dada por (5.13) y como 
función de rigidez sólida efectiva, se utiliza la función a' e (u) dada por (5.14), donde los 
parámetros necesarios 

Voo = -2,7 x 10~ 4 [rjM -1 ], C = 21,5, u max = 0,5, a = 5,7[Pa] y n = 5, (5.38) 

corresponden al modelo de suspensión con compresión tipo Kaolin (ver [9]). Además Ag = 
1690 [Kg/m 3 ] y g = 9,81 [Kgm/s 2 ]. 

La figura 5.11 muestra las funciones modelo f{u) y a(u) para este caso. 




0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 



Figura 5.11: Funciones modelo f(u) (izquierda) y a(u) (derecha) para el problema de sedimenta- 
ción-consolidación, segundo ejemplo. Las unidades son [m/s] para f(u) y [m 2 /s] para a(u). 



En [9] se da la siguiente expresión para el término difusivo integrado: 

( 



A(u) 



0. 



SÍ U < U r 



A{u) — A(u c ), SÍ U > U c 



(5.39) 



donde 



Aya) = 1 



C n 



j=l \l=i 



K n + 1 — l\ fu m í 



C + l 



u 



1 , 



(5.40) 
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cuya gráfica se muestra a continuación. 




O 0.1 0.2 0.3 0.4 0.5 

11 

Figura 5.12: Término difusivo integrado A(u) para el problema de sedimentación-consolidación, 
segundo ejemplo. 

En la tabla 5.3 se muestran la proporción V, tasa de compresión y errores entre la 
solución obtenida utilizando multiresolución y la solución obtenida sin multiresolución. 



t[a] 


V 


A 4 


ei 




&2 








60 


1.4109 


5.6100 


4.31x10" 


5 


2.34x10- 


4 


1.46x10- 


4 


2000 (*) 


4.4782 


7.1542 


6.87x10- 


-5 


5.78x10- 


-4 


7.88x10- 


-4 


6000 (*) 


7.2384 


10.7245 


1.36x10- 


4 


9.45x10- 


4 


9.65x10- 


4 


10000 (*) 


10.4568 


10.9781 


6.74x10" 


4 


1.32x10- 


3 


1.03x10- 


3 



Cuadro 5.3: Caso batch de suspensiones floculadas, segundo ejemplo. Tolerancia prescrita e — 
10 -3 , N = 129 puntos en la malla fina y L = 5 niveles de multiresolución. (*): figuras 5.13 - 5.15. 

Análogamente al primer ejemplo, en la tabla 5.3 puede verse que los errores entre la 
solución obtenida utilizando multiresolución y la solución obtenida sin multiresolución, 
están por debajo de la tolerancia prescrita. De igual modo, se ve una gran rebaja en costo 
computacional, dada por la alta tasa de compresión y proporción V. 

En la figura 5.13 se presenta un perfil de concentración en í = 2000 [s], para la solución 
utilzando multiresolución, y la solución sin multiresolución. La solución se calcula utili- 
zando 129 puntos en la malla fina. Se presenta además la configuración correspondiente 
de los coeficientes de ondelette significativos. 
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Segundo ejemplo: sedimentación batclr eon compresión í =2000[s] Coeficientes de ondelcttc significativos , jV =129, L = 5, £ = 10^ 



0.16 
0.14 
0.12 
0.1 
x 0.08 
0.06 
0.04 
0.02 




* i=2000[s] (MR) 

— í=2000[s] (DF) 

- - u =110 
u—u c 




0.16 
0.14 
0.12 
0.1 
x 0.08 
0.06 
0.04 
0.02 
















+ 
+ 


+ 




+ 
+ 


+ 

+ 
+ 


+ 

+ 




+ 
+ 
+ 
+ 
+ 
+ 








+ 
+ 


+ 


+ 



1 



2 3 4 

Niveles de Multiresolucion 



Figura 5.13: Izquierda: Condición inicial (rayas) y perfil de concentración a í = 2000[s] para 
el problema de sedimentación-consolidación, segundo caso (asteriscos). Derecha: Coeficientes de 
ondelette significativos correspondientes, e = 10~ 3 , -/V = 129 y L = 5. 



Segundo ejemplo: sedimentación batch con compresión f =fi000[s] 
0.1 6 r 

* i =6000[s] (MR) 

1 =6000[s] (DF) 

— — — u =TlQ 
u =u c 




0.05 0.1 

u(a;) 



0.15 



Coeficientes de ondelette significativos. N Q =129, L = 5, £ = 10 
0.16 




2 3 

Niveles de Mrdtiresolucion 



Figura 5.14: Izquierda: Condición inicial (rayas) y perfil de concentración a t = 6000 [s] para 
el problema de sedimentación-consolidación, segundo caso (asteriscos). Derecha: Coeficientes de 
ondelette significativos correspondientes, e = 10~ 3 , A^o = 129 y L = 5. 



En las figuras 5.14 y 5.15 se presentan perfiles de concentración en tiempos t — 6000[s] y t = 
10000[s], utilizando multiresolución, y la configuración de coeficientes de ondelette significativos. 
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Segundo ejemplo: sedimentación bateh con compresión t = 10000[s] Coeficientes de ondelette significativos . Nq = 129, L = b, £ = 10^ 

0.1 6 i — t ' ' ' — i — ' ' ' ' 1 0.1 6 1 1 1 1 1 1 1 1 1 

>\ | * t=10000[s] (MR) 




u(x) Niveles de Multare solución 



Figura 5.15: Izquierda: Condición inicial (rayas) y perfil de concentración a í = 10000[s] para 
el problema de sedimentación-consolidación, segundo caso (asteriscos). Derecha: Coeficientes de 
ondelette significativos correspondientes, e = 1CP 3 , -/Vq = 129 y L = 5. 



Finalmente se presenta en la figura 5.16 la solución numérica del problema de sedimentación 
consolidación obtenida utilizando el método de multircsolución, hasta el tiempo t — I2[h\. 



Perfiles de concentración hasta t — 12[h] 



x[m] 0.08 














































0.05 0.1 0.15 



Figura 5.16: Perfiles de concentración hasta t = 12[h] para el segundo problema de sedimentación- 
consolidación, asentamiento Batch. e = 1CP 3 , TVo = 129 y L = 5. 



Los resultados numéricos concuerdan con los resultados obtenidos por Bürger y Karlsen [9]. 
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5.3.4. Simulación de sedimentación continua 

Se modela un ICT de longitud 2, con una concentración inicial, uq = 0,052. En x = 1 
se prescribe una alimentación dada por *(í) = -8,55 x 10~ 7 . Se supone el ICT cerrado, 
es decir, q = y se simula el proceso de llenado hasta antes que el nivel de concentración 
en x = alcance el valor u(0, t) = 0,171. En ese momento, el recipiente se abre, y se hace 
q(t) = — 5 x 10 _6 [m/s]. Notar que desde ese momento, *(£) = 0,171 -q(t), es decir, el flujo 
en la alimentación es igual al flujo de descarga y el perfil de concentración entra en estado 
constante [9]. 

Notar que en este caso se utiliza como modelo el problema B (5.7)-(5.10). Se utiliza 
una función de densidad de flujo dada por 



f(u) = -1,98 x 10~ 4 u f 1 



0,3 



5,647 



(5.41) 



y una función de rigidez sólida efectiva dada por 



a e (u) 



0, 

5,7 



1 



si u ^ u c := 0,1, 
si u > u c := 0,1. 



(5.42) 



Estas aproximan a las funciones modelo determinadas para suspensión de carbonato de 
calcio [9]. 




3.5 
3 
2.5 

1.5 
1 

0.5 





Figura 5.17: Funciones modelo f(u) (izquierda) y a(u) (derecha) para para la simulación de 
sedimentación continua. Las unidades son [m/s] para f(u) y [m 2 /s] para a(u). 
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En este caso, Ag = 1690 [Kgm 3 ]. Además 
0. 

6,1267 x 10 2 • u 8 ( 1 - ¿ 



a(u) 



5,647 



si u sí Vb C := 0,1, 



si u > u c := 0,1, 



(5.43) 



y para el término difusivo integrado A(u) se utiliza (5.39), (5.40). Su gráfica se muestra 
en la figura 5.18. 



A(u) 




Figura 5.18: Término difusivo integrado A(u) para el problema de sedimentación continua. 



En la tabla 5.4 se muestran la proporción V, tasa de compresión y errores entre la 
solución obtenida utilizando multiresolución y la solución obtenida sin multiresolución. 



t[s] 


V 




ei 








eoo 




1800 


6.6818 


16.0156 


7.81x10- 


-5 


5.83x10- 


-5 


1.80x10- 


-5 


3600 (*) 


7.0845 


16.0156 


1.61x10" 


-4 


6.77x10- 


-5 


4.01x10- 


-5 


7200 (*) 


7.6731 


15.3010 


2.44x10" 


-4 


9.05x10- 


-5 


6.46x10- 


-5 


14400 (*) 


9.5790 


14.6441 


4.92x10- 


-4 


1.64x10- 


-4 


1.84x10- 


-4 


43200 (*) 


14.0489 


19.6441 


5.10x10- 


-4 


4.26x10- 


-4 


4.76x10- 


-4 



Cuadro 5.4: Simulación de sedimentación continua. Tolerancia prescrita e = 5 x 10 4 , A^o = 513 
puntos en la malla fina y L = 5 niveles de multiresolución. (*): figuras 5.19 - 5.22. 

Al mirar la tabla 5.4, de nuevo los errores entre la solución obtenida utilizando mul- 
tiresolución y la solución obtenida sin multiresolución, se encuentran por debajo de la 
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tolerancia prescrita. Una alta tasa de compresión y proporción V de tiempo total de CPU 
delatan la importancia del método de multiresolución en la aplicación de este tipo de 
problemas. 
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Figura 5.19: Izquierda: Condición inicial (rayas) y perfil de concentración a, t — l[h) para el 
problema de sedimentación continua (asteriscos). Derecha: Coeficientes de ondelette significativos 
correspondientes. e = 5x 10 -4 , Nq — 513 y L = 5. 



Simulación de sedimentación continua í=7200[s] Coeficientes de ondelette significativas, Ai, =513, i = 5, e=5 xlO" 3 
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Figura 5.20: Izquierda: Condición inicial (rayas) y perfil de concentración ai = 2[h] para el 
problema de sedimentación continua (asteriscos). Derecha: Coeficientes de ondelette significativos 
correspondientes, e = 5 x 1CP 4 , Nq = 513 y L = 5. 



En la figura 5.22 se presenta un perfil de concentración para el modelo de sedimentación 
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Simulación de sedimentación continua t = 14400[s] 
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Figura 5.21: Izquierda: Condición inicial (rayas) y perfil de concentración a t = 4[/i] para el 
problema de sedimentación continua (asteriscos). Derecha: Coeficientes de ondelette significativos 
correspondientes, e = 5 x 1CP 4 , iY = 513 y L = 5. 



continua, a í = 43200[s]. Notar que en este tiempo la solución ya entra en un estado 
estacionario, pues el flujo de alimentación es igual al flujo de descarga. 
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Figura 5.22: Izquierda: Condición inicial (rayas) y perfil de concentración ai— Í2[h] para el 
problema de sedimentación continua (asteriscos). Derecha: Coeficientes de ondelette significativos 
correspondientes, e = 5 x 1CP 4 , N = 513 y L = 5. 



76 



5.3 Ejemplos numéricos 



77 



Finalmente se presenta en la figura 5.23 la solución numérica del problema de sedi- 
mentación continua, obtenida utilizando el método de multiresolución, hasta el tiempo 
í = 16[h]. 



Distintos tiempos de simulación de sedimentación continua, hasta t = 16[/¿] 




0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 



Figura 5.23: Perfiles de concentración hasta t = 16[h] para el problema de sedimentación continua. 
e=5x 1(T 4 , N = 513 yL = 5. 



Los resultados numéricos concuerdan con los resultados obtenidos por Bürger y Karlsen 



[9]. 
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Capítulo 6 



Conclusiones y perspectivas 



6.1. Conclusiones 

En el presente trabajo se desarrolló un esquema numérico completamente adaptativo 
para acelerar los cálculos de volúmenes finitos de ecuaciones diferenciales parabólicas (ori- 
ginalmente desarrollado para leyes de conservación hiperbólicas) y ecuaciones parabólicas 
fuertemente degeneradas en una dimensión espacial. Se estudiaron varios casos test de 
ecuaciones hiperbólicas, parabólicas finales y no lineales, y ecuaciones parabólicas fuerte- 
mente degeneradas provenientes de la teoría de procesos de sedimentación-consolidación. 

Generalmente, al añadir un término viscoso a un esquema, la solución tiende a suavizar 
y en algunos casos puede estabilizar un esquema numérico originalmente inestable. Se pudo 
ver que excepto por una limitación de paso temporal (que en el caso invíscido es diferente) 
el problema viscoso no implica mayores complicaciones desde el punto de vista numérico. 

El análisis de multiresolución se mantiene inalterado, pues sólo tiene que ver con la 
regularidad de los valores puntuales o medias en celda de la solución. 

Es importante destacar que en el capítulo 5 se utilizaron esquemas de diferencias finitas, 
por lo que en los algoritmos de multiresolución empleados se considera un análisis de 
multiresolución para valores puntuales. 

Se comienza con una discretización de volúmenes finitos (o diferencias finitas) en una 
malla uniforme, y una integración explícita en tiempo, ambas de segundo orden. Mediante 
técnicas de análisis de multiresolución, se reduce el tamaño de la malla, eliminando los 
puntos con detalles no significativos, pero manteniendo siempre un esquema de segundo 
orden. 

La actualización temporal de la malla se realiza mediante una estrategia de adaptación 
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dinámica que aprovecha la representación puntual esparsa, agregando coeficientes vecinos 
en escala y espacio para mejorar la captura de la información. 

Para la evaluación de los flujos numéricos, en la malla localmente refinada, se utilizaron 
esquemas ENO de segundo orden y esquemas de Engquist-Osher modificados de segundo 
orden. 

Los algoritmos empleados son generalizables al caso de otras condiciones de borde 
(simplemente modificando el interpolador intermallas y el cálculo de los flujos en los puntos 
de frontera), otra elección para la condición inicial, otro tipo de evolución temporal, otra 
elección para los predictores intermallas, otra elección para el orden de las interpolaciones 
ENO, otra elección para el cálculo del flujo numérico, otro tipo de estructura de datos, 
etc. 

La eficiencia del algoritmo fue medida mediante la tasa de compresión y el tiempo de 
CPU. La diferencia de tiempo total de CPU entre la solución numérica que no utiliza 
multiresolución y la que utiliza multiresolución está directamente relacionada con el hecho 
de que en una, la solución numérica sin multiresolución se evalúan todos los flujos numéri- 
cos mientras que en la otra solución numérica con multiresolución, sólo se calculan los 
flujos numéricos donde existen coeficientes de ondelette significativos. Lógicamente esta 
diferencia se ve incrementada cuando el flujo numérico es más costoso. 

La aplicación del método de multiresolución resulta aún más provechosa en la simula- 
ción de procesos de sedimentación de suspensiones floculadas. El que las ecuaciones sean 
de naturaleza más compleja, se suma el hecho de que los resultados experimentales publi- 
cados requieren un tiempo de simulación de varias horas, en contraste con las fracciones 
de segundo suficientes para estudiar la solución numérica de los problemas hiperbólicos y 
parabólicos incluidos en este trabajo. Además, la condición CFL en este caso, hace que 
Ai sea muy pequeño. Esto hace pensar en la utilización de un esquema implícito o semi- 
implícito [6]. 

La gran desventaja de utilizar algoritmos de multiresolución, es quizás el hecho de que 
los resultados en cuanto a convergencia aún no tienen un gran auge. Una gran parte de 
los argumentos del análisis de multiresolución desarrollado por Harten es de naturaleza 
heurística. 

En la parte final se presentó un método numérico para obtener soluciones aproximadas 
de problemas provenientes de fenómenos de sedimentación. La idea desarrollada fue aplicar 
los métodos de multiresolución a los esquemas diseñados por Bürger et al. [5, 7, 8, 9, 10] 
y se observó que el método de multiresolución es de gran ayuda para reducir el costo 
computacional en este tipo de problemas sin afectar la calidad de la solución. 

Todos los experimentos se realizaron en equipos con procesadores Pentium 4 de 1.6 
Mhz, con 1GB de memoria RAM, tanto en plataforma Linux como Windows. 
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6.2. Perspectivas 

■ Para el caso de ecuaciones parabólicas fuertemente degeneradas, la perspectiva a 
más corto plazo es modificar el algoritmo para poder aplicarlo a las ecuaciones que 
modelan otros tipos de fenómenos de sedimentación. 

■ Aplicar métodos de multiresolución a la resolución de problemas inversos. 

■ Utilizar esquemas ENO de orden superior a dos. Combinar esto con la utilización de 
esquemas con varios switches [24]. 

■ Aplicar métodos de multiresolución a problemas que modelan la separación de sus- 
pensiones polidispersas [4]. 

■ Realizar los experimentos del capítulo 5 utilizando esquemas semi-implícitos e implíci- 
tos. Esto se traduce en pasar de un Ai de orden de (Ax) 2 a un orden de Ax. Sin 
embargo las complicaciones están en tener que resolver un sistema de ecuaciones 
no-lineales en cada iteración. Además el proceso de multiresolución para esquemas 
semi-implícitos se complica bastante. 

■ Extender los resultados de los puntos anteriores al caso de sistemas y ecuaciones 
multidimensionales . 

■ Los códigos pueden ser fácilmente traducidos a un lenguaje más robusto como 
FORTRAN, C, o C+- 1-, dado que las funciones y subrutinas en la implement ación no 
abusan de las funciones implícitas de MATLAB (excepto en la estructura SPARSE 
de los datos). 
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Apéndice A 



Cálculo de los coeficientes de 
interpolación en la multiresolución 

A.l. Multiresolución de valores puntuales 

En la sección correspondiente se ha mencionado que 

s 

l{x k 2 j\, U k ) = A(«J+í-i + uj-i) (A.l) 
1=1 

es el polinomio de grado r — 1 que interpola los puntos (v,j_ s , . . . , u k +s _ l ). Para ver esto, y 
encontrar los valores de los coeficientes (3i , se utiliza el polinomio interpolador de Lagrange 

j+s-l j+s-l k 

l=j—s t=j—s 1 1 

en el punto £ 2 j-i 

i(x k j} v u k )= e «? n ^ 

l=j—s l=j—s l 1 

donde xí^-i = (2j — 1) ■ = (j -1/2) -h k . Luego 

j+s-l j+s-l ■ _ 1 _ ■ 

I 7 tí 

(=j-s <=j-s 
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Si se toma en cuenta que los pares de valores puntuales (uj-i, Uj), (uj^Uj+i), • • • están 
multiplicados por el mismo factor, se tiene (A.l), con 

_ l 2 -3 2 ---(2¿-3) 2 -(2¿-l)-(2¿ + l) 2 ---(2. S -l) 2 m 
Pl 2 2 *" 1 ■ (s + l- 1)1 ■ (s-l)l { ' ' 

Por lo tanto se tienen los siguientes coeficientes para cada r = 2s mencionado: 
- r = 2, 8 = 1, = 2ok • (-1) 2 = §• 

■ r = 4, s = 2, = ¿fjy • (-1) 2 = ±fo= X¿ {)[ ■ (-1) 3 = 
A. 2. Multiresolución de medias en celda 

De manera análoga al caso anterior, en la sección correspondiente se ha mencionado 
que 

d k- ü k-i _ ü k-i _ ü k-l I ( x 2j-v uk ) ~ U j-i (A2) 

a j- U 2j-l U 2j~l~ U 2j~l . i \ K - ¿ ) 



Se aplica el caso anterior (para valores puntuales) al esténcil (Uj_ 8 , . . . , Uj +S _ 1 ), por tanto 

s 

1{x\~\- u k ) = Y,Wt+i-i + U^i), 
1=1 

con los mismos f3¡ calculados en el apéndice A.l, por tanto 

,s fí (TTk i Tfk \ Jjk 



, x _ELiA(^-i + ^)-^-i 
Ur >,- , 

i 



h k . 



y utilizando hk = 2h k -i, se tiene 



u 



fc-l 



2MU[ + Uj-i) + 202(U} +1 + Cf*_ 2 ) + ■ ■ ■ + 2P S (U^ +S _ 1 + U}_ a ) - 2\J]_ X 

hk 

o equivalentemente, 

h k • ñ^_\ = • • • + (2A - 1 + 2(3 2 ) ■ (U}_ 2 - U}_ 3 ) + (2/3x - 1) • (U^ - U}_ 2 ) + 
1 • (U* - U}_{) + (1 - 2/30 • " C?) + 

(l-2/3 1 -2/3 2 )-(C/jV 2 -^i) + --- 
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Además, si se toma en cuenta la relación 



j ~ h k ' 



y que las medias en celda u k - +l y Uj +l , l = 1, . . . , s — 1, están multiplicados por el mismo 
factor (sólo cambia de signo), se llega a la expresión 



s-1 

-k 



í=i 

con 

ji = -(2-Pi -ji-i), 70 = 1- 
Por lo tanto se obtiene para cada f = 2s — 1 mencionado: 



7i = -(2-/3 1 -7ü) = -(2-^-1) = -|. 



r = 3, s = 2, 

_ — _fo . _ 

16 

f = 5, s = 3, 

71 = _(2./3 1 - 70 ) = -(2.1i-l) = -^, 

72 = -(2-/3 2 -7i) = -(2-if + ^ = 4. 
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Apéndice B 



Análisis de estabilidad para el caso 
parabólico no lineal 



A continuación se analizará la estabilidad en el sentido de la variación total del esquema 
numérico presentado en la sección 4.1.7. Este análisis puede aplicarse al caso de ñujos 
lineales o no lineales. Se quiere encontrar una condición CFL que pueda utilizarse para el 
esquema ENO-TVD de segundo orden. 

Un esquema se dice TV-estable si la variación total 



N-l 



TV(v h (x,t)) = TV(v n ) := ^ \v] +l - v?\ 

3=0 

de una sucesión de aproximaciones numéricas Vh(x,t) está acotada uniformemente en 
h = Ax y t = nAt, con /i^OyO^i^T. Aún más, el esquema es TVD si 

TV{v n+l ) ^ TV{v n ). 

Claramente, un esquema TVD es TV-estable. 

Con estas definiciones básicas, Harten [25] probó el siguiente 

Lema 1 Si un esquema escrito en la forma 



v 



n+1 

3 



+ C1-A+V? -Cj^A^v], (B.l) 



satisface, para todo j, 



C+ > 0, (B.2) 
CJ > 0, (B.3) 
C+ + CT ^ 1, (B.4) 
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entonces el esquema es TVD. 

Él introduce el siguiente esquema explícito, de primer orden, 

(B.5) 
(B.6) 

(B.7) 
(B.8) 

donde fj = f(uj), y gj = g(uj) es elegida tal que 

I&KpK-) (b.9) 

donde p es el clásico limitador de flujo de Harten 



,n+l 

i 














= ^ 
A-Vj 






A +9j 





p(a) 



0, para primer orden en espacio 

|(|o| — a 2 ), para segundo orden en espacio. 



Con estas definiciones, Harten prueba que para esquemas de primer y segundo orden, una 
condición suficiente para que el esquema sea TVD es la condición tipo CFL 

máx|tjj| 1, (B.10) 

3 

pues ujj es el coeficiente CFL medio local. 

Se quiere modificar la demostración hecha por Harten [25] para el caso de esquemas 
de segundo orden, con el fin de aplicarla al caso viscoso, para ello, Bihari [3] probó el 
siguiente 

Teorema 3 Un esquema escrito en la forma (B.l), con C^ 1 definido por 

c t = ó t + x h (B11) 

es TVD si 

a < — — (B.12) 
Re + 4 v ' 

con 

Ax 

a = máx \u)j\, Re = máx\ujj\— — . (B.13) 
j 3 
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Notar que el esquema (B.l), (B.ll) es de segundo orden en espacio al aproximar la solución 
del problema, puesto que se ha incluido un término viscoso con una discretización central 
al esquema original TVD de segundo orden. 

Notar además que la definición dada para a y Re difieren de las definiciones dadas en 
(4.24) y (4.25) para el caso lineal. Sin embargo el significado cualitativo de estas cantidades 
es el mismo. Es decir, (B.13) es la definición equivalente para a y Re en el caso no lineal. 

Dem: (Del teorema) Se mostrará que se satisfacen las condiciones del Lema 1. Con la 
definición dada de (B.6), es claro que se satisfacen las condiciones (B.2) y (B.3). Falta 
entonces mostrar que 



De las definiciones dadas y de la propiedad (B.9) se sigue que (B.14) se satisfará si 




(B.14) 




Ahora, dado que a ^ 1 (necesario para que se satisfaga (B.9)), se tiene que a^(a-l) ^ 0. 
Luego, es posible obtener una versión levemente más restrictiva que (B.15): 




la cual se satisface si se satisface 



Re 



cr s$ 



Re + 4 



□ 
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Apéndice C 



Código y documentación 



Tanto los códigos en MATLAB para cada experimento, la documentación respectiva, 
como una versión electrónica de este informe pueden ser obtenidos en forma gratuita, 
desde el sitio http:/ /www.udec.cl/~riruiz/tesis.html. 
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